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ABSTRACT 


This  dissertation  provides  experimental  observations  and  analyses  that 
associate  low-latitude  transionospheric  signal  scintillation  with  transequatorial  VHF 
radio  propagation  and  errors  in  transionospheric  geopositioning. 

The  experiment  observed  equatorial-region  ionospheric  total  electron  content 
(TEC)  derived  from  Global  Positioning  System  (GPS)  signals  using  receivers  on 
Oahu,  Hawaii,  Christmas  Island,  and  Rarotonga,  Cook  Islands.  The  experiment 
simultaneously  measured  VHF  transequatorial  propagation  of  VHF  television  signals 
from  Hawaii  to  Rarotonga. 

Analysis  shows  that  a  moving  second  moment  of  vertical-equivalent  TEC 
strongly  correlates  to  each  VHF  transequatoral  radio  propagation  event.  From 
experimental  observation  analysis,  the  author  develops  models  for  prediction  of  TEP 
and  time-space  distribution  of  low-latitude  transionospheric  scintillation. 

The  author  also  develops  equations  that  show  the  potential  errors  in  time, 
frequency,  and  angle  used  in  geopositioning  solutions.  These  three  parameters  are 
potentially  correctable  using  these  techniques. 
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I.  INTRODUCTION 

The  reflecting  layers  of  the  ionosphere  have  been  used  to  support  long- 
distance communications  since  the  early  days  of  radio.  In  past  years,  many 
communicators  have  used  HF  (2  to  30  MHz)  for  long  distance  communications,  and 
this  band  is  still  widely  used  for  this  purpose.  Predictions  of  HF-propagation 
conditions  have  evolved  to  help  communicators  contend  with  diurnal,  seasonal,  and 
sunspot-cycle  variations  of  the  ionosphere.  Recent  advances  in  microcomputers  have 
allowed  HF-propagation  prediction  programs  to  become  available  to  a  large  number 
of  users  at  a  very  low  cost. 

While  HF-prediction  programs  are  highly  useful  to  mid-latitude  users,  they 
have  many  limitations.  At  high  latitudes  (above  about  60°)  and  low  latitudes  (below 
about  25°),  prediction  programs  are  less  reliable.  The  mid-latitude  ionospheric 
models  do  not  adequately  describe  the  ionospheric  behavior  of  these  areas  (the  two 
polar  regions  considered  to  be  one  genre).  These  areas  exhibit  anomalous 
propagation  and,  as  such,  are  unpredictable.  However,  through  anecdotal  reports 
from  amateur  radio  operators  and  other  observers  of  radio  propagation,  it  has  been 
noted  that  these  areas  not  only  support  HF  communications  but  also  VHF 
communications.  Scientific  exploration  has  contributed  much  to  the  knowledge  of 
the  morphologies  of  these  regions,  but  attempts  at  communications-path  predictions 
have  generally  been  unsuccessful.  Successful  prediction  of  nighttime  transequatorial 
VHF  communications  is  the  first  of  two  major  objectives  of  this  dissertation. 


Radio  signal  paths  which  originate  in  the  northern  and  terminate  in  the 
southern  hemisphere  at  a  near-conjugate  point  (and  vice  versa)  are  known 
collectively  as  transequatorial  propagation  (TEP)  paths.  TEP  is  further  divided  into 
two  types,  identified  as  afternoon-TEP  and  evening-TEP.  Afternoon-TEP  is 
characterized  by  chordal  propagation  from  one  equatorial  anomaly  to  the  other  and 
thence  return  to  earth.  Afternoon-TEP  is  generally  time-restricted  to  between  about 
1000  local  and  sunset  and  supports  frequencies  up  to  about  60  MHz.  Similar  to  mid- 
latitude  propagation,  the  probability  of  afternoon-TEP  increases  with  increasing 
sunspot  numbers.  Evening-TEP  begins  sometime  after  sunset  when  signals  propagate 
via  an  ionospheric  ducting  mechanism.  Since  the  signals  are  propagated  within 
ducts,  frequencies  and  strengths  of  the  signals  are  higher  than  would  otherwise  be 
expected.  Reports  of  evening-TEP  propagation  at  frequencies  as  high  as  432  MHz 
have  been  reported  (Cracknell  and  Whiting,  1980),  and  reception  at  frequencies  as 
high  as  200  MHz  was  detected  in  this  research  effort. 

The  effects  on  signals  influenced  by  the  ducting  mechanism  mentioned  are  not 
all  enhancements,  however.  Severe  degradation  of  signals  which  pass  through  the 
ducts,  i.e.,  to  and  from  satellites,  may  prevent  their  reception  and  demodulation 
altogether.  The  signal  degradation  is  a  result  of  amplitude,  phase,  and  frequency 
scintillation  ±at  is  induced  as  the  signal  passes  through  the  duct.  One  process  which 
utilizes  these  transionospheric  signals  for  its  operation  and  is  significantly  degraded 
by  the  ducts,  is  the  geolocation  of  radio  transmitters. 


Geolocation  can  be  accomplished  by  making  measurements  at  the  receiver  site 
of  angle  of  arrival  (AOA),  time  difference  of  arrival  (TDOA),  and  frequency 
difference  of  arrival  (FDOA).  Both  AOA  and  TDOA  systems  locate  the  emitter  using 
multiple  bearing  determinations  while  FDOA  systems  use  multiple  range 
determinations. 

If  the  receiver  is  aboard  a  space-borne  vehicle,  and  the  emitter  is  on  or  near 
the  earth,  the  electromagnetic  signal  must  pass  through  the  ionosphere  on  its  journey 
from  the  emitter  to  the  receiver.  Errors  are  introduced  by  the  ionosphere  to  the 
measurement  of  angle,  time,  or  frequency  of  arrival  at  the  receiver. 

If  the  ionospheric  region  through  which  the  signal  propagates  is  in  the  mid- 
latitudes  of  the  earth,  the  structure  and  density  of  the  ionospheric  layers  are 
sufficiently  known  and  predictable.  The  mid-latitude  ionosphere  (in  an  undisturbed 
condition)  can  therefore  be  modeled  with  enough  accuracy  to  compensate  for  the 
angle-  and  time-of-arrival  errors. 

In  both  high  and  low  latitudes  the  ionosphere  is  not  stable  nor  predictable. 
At  present,  no  reliable  model  exists  to  predict  ionospheric  behavior  in  these  regions. 
With  no  reliable  models,  the  errors  introduced  to  the  geolocation  factors  cannot  be 
compensated  for,  and  the  location  of  the  emitter-of-interest  may  be  unobtainable. 
Characterization  of  scintillation  induced  by  the  low-latitude  ionosphere  is  the  second 
major  objective  of  this  study. 

The  ionosphere  of  the  low-latitudes  is  introduced  in  Chapter  II.  The 
equatorial  E-layer  electric  and  magnetic  fields  are  described,  which  combine  to 
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produce  the  equatorial  anomalies.  The  instability  of  the  dissipating  anomalies  is 
recounted  along  with  duct  formation.  Methods  for  observing  the  ionosphere 
remotely  are  detailed  in  Chapter  III  with  emphasis  on  equipment  which  was  used  in 
this  research  campaign.  The  theories  of  low-latitude  transequatorial  and 
transionospheric  radiowave  propagation  are  developed  in  Chapter  IV.  Chapter  V 
contains  the  experimental  setup  whereby  two  oppositely-directed  radio  transmitter- 
receiver  links  were  established  in  the  mid-Pacific  with  accompanying  ionospheric 
monitors.  And  Chapters  VI,  VII,  and  VIII  include  analysis  of  the  collected  data  and 
derivations  of  the  prediction  algorithms.  Algorithms  developed  in  these  chapters  are 
included  in  the  appendices  as  MATLAB  programs. 

The  material  of  the  dissertation  falls  into  two  major  segments.  The  initial, 
Chapters  I  through  IV,  is  considered  background  and  explanatory.  The  remaining, 
given  in  Chapters  V  through  IX,  contains  original  and  novel  work  including: 

*  Derivation  of  low-latitude  ionospheric  predictors  based  upon  mid- 
Pacific  experiment 

*  Prediction  of  TEP  based  upon  TEP  observations 

*  Correlation  of  TEP  and  GPS-derived  total  electron  content 

*  Low  latitude  scintillation  effects  upon  geolocation 

*  Temporal  and  spatial  scintillation  models  using  GPS-derived  total 
electron  content,  and 

*  Bi-directionality  of  TEP  circuit. 


II.  THE  PROPAGATION  ENVIRONMENT 

Free  ions  and  electrons  existing  above  about  90  km  altitude  constitute  the 
earth's  ionosphere.  At  low  latitudes,  on  each  side  of  the  magnetic  equator,  the 
structure  of  the  ionosphere  differs  from  any  other  location.  The  magnetic  equator 
is  the  only  location  in  the  terrestrial  ionosphere  where  both  the  electric  and  magnetic 
fields  are  parallel  to  the  surface  of  the  earth.  The  effects  of  this  singularity  create  the 
low-latitude  ionosphere  with  its  equatorial  anomalies.  The  cause  of  this  deviation 
is  the  complex  interaction  between  the  earth's  magnetic  field  and  the  motions  of 
electrons  and  ions  in  the  equatorial  ionospheric  region.  Understanding  of  the 
physics  of  the  environment  is  the  objective  of  this  chapter. 
A.         SOLAR  EMISSIONS 

Emissions  from  the  sun  provide  the  energy  necessary  to  ionize  the  neutral 
molecules  of  the  upper  atmosphere.  The  ionizing  energy  arrives  in  two  forms: 
electromagnetic  (as  either  x-rays  or  ultraviolet  rays)  and  charged  particles  (electrons, 
protons,  and  ions). 

Radio  noise  is  also  emitted  from  the  sun.  The  power  density  of  noise  at  2800 
MHz  (wavelength  10.7  cm)  is  measured  daily  at  a  receiving  site  near  Ottawa.  This 
power  density,  called  the  10.7  cm  flux,  given  in  units  of  10"22  W  m 2  Hz"1,  is  reported 
to  scientists  at  the  National  Geophysical  Data  Center  (NGDC)  in  Boulder,  CO.  They 
in  turn  report  daily  readings  (and  past  values)  on  their  Internet  homepage  whose 


address  is  http://www.sel.bldrdoc.gov/today.html.  The  flux  values  vary  from  about 
65  to  250  units  at  the  nadir  and  peak  of  the  solar  cycle,  respectively. 

Solar  emissions  are  directly  proportional  to  the  sunspot  number.  This  is  due 
to  heightened  flare  activity  connected  with  an  increase  in  the  number  of  sunspots. 
Flares  emit  electromagnetic  energy  creating  more  x-rays  and  ultraviolet  rays. 

The  10.7  cm  solar  flux  also  increases  with  the  elevated  flare  activity  and  as 
such  is  a  dependable  predictor  of  sunspot  number.    It  reliably  indicates  energy 
available  for  ionization  in  the  atmosphere. 
B.         THE  EARTHS  MAGNETIC  FIELD 

The  magnetic  field  of  the  earth  can  be  approximated  to  first  order  as  a  dipole 
field,  with  the  poles  tilted  about  11°  from  the  axis  formed  by  the  geographic  poles. 
The  north  geomagnetic  pole  is  located  approximately  at  geographic  latitude  and 
longitude  79.19°  N  and  70.98°  W,  with  the  south  geomagnetic  pole  located  at 
79.19°  S  and  109.02°  E.  The  magnetic  lines  of  flux  leave  the  south  magnetic  pole 
and  stream  toward  the  north  magnetic  pole. 

The  inclined  poles  create  a  need  for  a  new  magnetic  coordinate  system, 
complete  with  equator,  latitudes,  and  longitudes.  Since  the  equatorial  ionosphere 
is  somewhat  symmetric  about  the  magnetic  equator,  it  is  necessary  to  convert  from 
geographic  to  the  geomagnetic  coordinates  for  terrestrial  points  of  interest.  A 
method  to  compute  the  geomagnetic  coordinates  for  a  point  of  interest  is  given  in 
Appendix  A. 


1.         Lines  of  Flux 

Magnetic  lines  of  flux  originate  from  the  north  pole  of  a  dipole  magnet  and 
terminate  on  the  south  pole  (indicating  that  the  earth's  magnetic  poles  are 
misnamed).  As  with  a  bar  magnet,  however,  other  lines  originate  from  different 
locations  toward  the  northern  end  and  terminate  at  conjugate  points  on  the  southern 
end.  This  is  depicted  for  the  earth  in  Figure  1,  where  it  is  seen  that  ±e  lines  of  flux 
are  parallel  to  the  earth's  surface  over  the  magnetic  equator.  The  distance  from  the 


Figure  1.  Magnetic  lines  of  Flux, 
earth's  center  to  a  point  on  a  magnetic  field  line  is  found  as  r  =  R^  L  cos20,  where 
R^  is  the  earth's  radius,  L  is  the  height  of  the  line  at  its  equatorial  plane  crossing 
measured  in  earth  radii,  and  O  is  the  magnetic  latitude  of  the  point  of  interest.  The 
length  of  a  field  line  between  a  point  in  one  hemisphere  and  its  magnetic  conjugate 
in  the  opposite  hemisphere  is  given  by  (Davies,  1990) 


CG  +  I  ln(C  +  G)  -  -  In  L 

2 


(1) 


where  C2  =  3(L  - 1)  and  G2  =  4L  -  3. 

2.        Magnetic  Declination  Angle 

The  angle  from  a  point  of  interest  to  the  geographic  north  pole  is  by  definition 
360  degrees  or  due  north,  which  is  called  true  north.  The  great  circle  heading  from 
this  same  point  of  interest  to  the  magnetic  north  pole  is  magnetic  north.  The 
difference  between  true  north  and  magnetic  north  is  the  magnetic  declination  angle. 
A  method  to  compute  the  magnetic  declination  angle  is  derived  in  Appendix  B. 
C.         THE  EQUATORIAL  IONOSPHERE 

The  earth's  ionosphere  is  created  by  energy  from  the  sun  being  imparted  to 
neutral  molecules  in  the  upper  atmosphere.  When  this  energy  is  absorbed  by  a 
molecule,  an  electron  is  separated  from  the  molecule  leaving  a  free  electron  and  an 
ion.  The  density  of  free  electrons  is  the  fundamental  quantity  from  which 
propagation  characteristics  are  derived.  Prediction  or  measurement  of  the  electron 
density  is  essential  to  the  quantitative  characterization  of  the  ionosphere. 

Electron  density  is  generally  directly  related  to  the  solar  zenith  angle,  %.  The 
solar  zenith  angle  is  the  angle  between  a  line  drawn  directly  upward  from  the  point 
of  interest  on  earth  to  the  line  drawn  from  that  same  point  to  the  sun,  as  shown  in 
Figure  2.    An  accurate  method  of  computing  the  solar  zenith  angle  is  given  by 
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Bourges  (1985).  Since  free  electron  production  is  maximum  when  the  sun  is  directly 
overhead,  the  densities  are  generally  highest  in  summer  at  noon. 
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Figure  2.  Solar  zenith  angle. 


This  generality  does  not  hold  at  low  latitudes.  Directly  over  the  equator  the 
electron  density  is  actually  less  than  the  densities  north  and  south  of  the  equator  due 
to  the  equatorial  anomalies.  These  equatorial  (or  Appleton)  anomalies  are  caused  by 
upward  electron  movement  above  the  equator  due  to  interaction  between  the  earth's 
equatorial  magnetic  and  electric  fields.  The  equatorial  electric  field  is  discussed  next. 

1.         The  Equatorial  Electrojet 

The  winds  of  the  upper  stratosphere  interact  with  the  plasma  of  the  lower 
ionosphere.  The  prevalent  winds  in  this  region  are  directed  upward  during  the  day 
and  downward  at  night  (Woodman,  1970).  Electrons  and  ions  created  from  solar 
energy  are  swept  up  by  these  winds  and  move  initially  in  the  direction  of  wind  flow, 
as  they  attempt  to  cross  lines  of  flux  of  the  earth's  magnetic  field.    A  force 
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perpendicular  to  the  wind  direction  is  applied  to  the  charged  particles  by  the  lines 
of  flux.  This  results  in  the  ions  being  deflected  in  one  direction  and  the  electrons  in 
the  other  (i.e.,  one  west  and  the  other  east).  The  deflecting  force  due  to  the 
magnetic  field  is 

F  =  qV  x  B,  (2) 

where  F  is  the  force,  q  is  the  amount  of  charge,  V  is  the  velocity  of  the  charged 

particle,  and  B  is  the  magnetic  flux  density.  The  net  effect  of  the  deflections  is  a 
current  flow  directed  in  the  east- west  direction.  This  current  flow  is  called  the 
equatorial  electrojet. 

The  electrojet  occurs  approximately  105  km  above  the  earth  in  the  E  region. 
The  current  density  at  midday  attains  a  value  of  about  14  A/km2  west-to-east  which 
creates  an  electric  field  directed  in  the  same  direction  as  the  current,  peaking  at 
about  0.5  mV/m  (Forbes,  1981;  Hargreaves,  1992).  Therefore,  the  earth's  magnetic 
and  electric  fields'  interaction  with  the  neutral  winds  and  the  ionized  particles  of  the 
lower  ionosphere  generates  an  electric  field  which  encircles  the  globe  above  the 
magnetic  equator.  The  field  is  directed  from  both  the  day  and  the  night  sides  toward 
the  dusk  terminator,  as  shown  in  Figure  3. 

2.         The  Equatorial  Anomalies 

Above  the  magnetic  equator  there  are  two  fields,  both  parallel  to  the  surface 
of  the  earth.  The  magnetic  field  is  aligned  magnetically  south  to  north  and,  at  all 
points  above  the  magnetic  equator,  parallel  to  the  earth's  surface.  During  the  day, 
the  electric  field  is  directed  from  magnetic  west  to  magnetic  east  along  the  magnetic 
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Figure  3.    The  earth's  equatorial 
electric  field  (After  Kelley,  1989). 

equator  and  at  all  points  parallel  to  the  earth's  surface.  The  interaction  between 

these  perpendicular  fields  and  the  free  electrons,  which  are  created  by  the  ionizing 

energy,  produces  a  condition  called  the  fountain  effect. 

The  fountain  effect  is  the  electrodynamic  lifting  of  free  electrons  above  the 

magnetic  equator.  The  underlying  physics  can  be  seen  by  first  dividing  both  sides 

of  Eq.  2  by  the  charge,  q,  to  obtain  the  electric  field  E,  i.e., 


E  =  -  =  -V  x  B, 


(3) 


where  the  minus  sign  indicates  the  velocity  of  the  electrons.   Equation  3  can  be 
crossed  into  B  to  obtain 


E  x  B  =  (-V    x  B)  x  B, 


(4) 


where  the  velocity  which  is  perpendicular  to  both  E  and  B  is  Vi .  Application  of  a 
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common  vector  identity  to  Eq.  4  produces 

E  X  B  =  V  (B  •  B)  -  B(B  •  VJ 

(5) 
=  VB2. 

Thus  the  electrons  move  in  the  direction  of  E  crossed  into  B,  with  velocity 

t7        E  x  B 

v,  =  —7-  (6) 

B2 

It  is  observed  that  the  perpendicular  electron  velocity  is  independent  of  both 
its  mass  and  charge.  Indeed,  a  parallel  derivation  of  the  perpendicular  velocity  for 
the  ions  results  in  an  identical  equation.  For  the  earth's  orientation  of  the  magnetic 
and  electric  fields,  the  perpendicular  velocity  is  vertical:  upward  during  daylight  and 
downward  during  darkness.  A  representative  magnitude  of  the  midday  velocity  can 
be  found  by  assuming  the  electric  field  to  be  0.5  mV/m  and  the  magnetic  field  at  50 
/xWb/m2  to  produce  a  typical  initial  velocity  of  10  m/s  directed  upward. 

During  the  day  electrons  continue  to  be  lifted  until  the  effects  of  pressure, 
gravity,  and  the  magnetic  field  overcome  the  effects  of  the  E  x  B  field.  The  net 
result  of  the  electric  and  magnetic  field  interaction  is  that  free  electrons  and  ions 
over  the  magnetic  equator  are  forced  straight  up,  rising  to  altitudes  of  up  to  1000 
km  (Hargreaves,  1992).  The  electrons  and  ions  then  follow  the  magnetic  lines  of 
flux  back  toward  the  earth,  increasing  the  electron  density  both  north  and  south  of 
the  magnetic  equator  at  latitudes  of  about  10  to  20  degrees.  The  electron  density 
over  the  equator  is  thereby  diminished  due  to  the  fountain  effect.  The  fountain  effect 
is  illustrated  in  Figure  4. 
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Figure  4.  The  fountain  effect. 


After  sunset  the  E-region  neutral  winds  shift  from  upward  to  downward.  With 
the  velocity  vector  downward,  the  electrojet  reverses,  flipping  the  polarity  of  the 
electric  field.  The  E  x  B  plasma  motion  described  by  Eq.  6  now  forces  the  electrons 
and  ions  downward,  pulling  plasma  from  the  overhead  area  directly  over  the 
magnetic  equator  and  evacuating  the  upper  regions.  The  high  plasma  density  of  the 
anomalies  begins  to  migrate  back  toward  the  magnetic  equator,  following  the 
magnetic  field  lines,  as  shown  in  Figure  5.  The  dynamics  of  this  plasma  migration 
are  unstable.  If  an  outside  perturbation  stimulates  this  instability,  the  ensuing 
upheaval  can  affect  the  nighttime  ionosphere  by  contributing  to  the  formation  of 
ionospheric  "bubbles"  or  "ducts"  also  known  as  Equatorial  Spread-F  (ESF). 

13 


1000 


200 


Magnetic 
Equator 

Figure  5.  Plasma  migration  from  anomalies  to  equator. 


3.        Equatorial  Spread-F 

Early  ionospheric  researchers  using  ionosondes  observed  the  occurrence  of 
a  nighttime  spread-F  on  equatorial  ionograms.  (Examples  of  ionograms  with  ESF  are 
shown  in  Chapter  III.)  Their  term  for  this  phenomenon,  Equatorial  Spread-F,  still 
survives  today.  Other  names  found  in  the  literature,  emphasizing  different  aspects 
of  ESF,  are  ducts,  bubbles,  and  plumes.  ESF  is  also  responsible  for  scintillations 
which  occur  on  satellite  communications  signals  at  low  latitudes.  Equatorial 
scintillations  have  been  observed  at  frequencies  as  high  as  7  GHz.  (Hargreaves, 
1992) 

The  causes  of  ESF  are  still  under  discussion  in  the  scientific  arena,  with  many 
theories,  both  linear  and  non-linear,  being  advanced  by  different  researchers.  It  is 
not  known  what  initiates  ESF,  nor  why  it  occurs  one  night  and  not  the  next.  It  is 
known  that  nighttime  long-range  transequatorial  propagation  of  VHF  signals,  the 
observation  of  equatorial  spread-F,  and  satellite  signal  scintillation  are  all  different 
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manifestations  of  the  same  physical  event.  The  onset  of  ESF  has  its  origin  in  the 
formation  of  ionospheric  ducts. 

4.         Equatorial  Ionospheric  Ducts 

After  sunset,  the  ionizing  energy  from  the  sun  vanishes  and  the  formation  of 
ions  and  free  electrons  in  the  ionosphere  ceases.  Equatorial  anomalies  lose  their 
source  of  electrons  and  begin  a  process  of  equalizing  the  densities  across  the  equator, 
as  was  shown  in  Figure  5.  Left  undisturbed,  anomaly  electrons  migrate  back  toward 
the  equator,  refill  the  pocket  of  diminished  electron  density  over  the  equator  and, 
sometime  after  local  midnight,  the  equatorial  ionosphere  is  indistinguishable  from 
the  mid-latitude  ionosphere.  This  is  precisely  what  occurs  on  some  evenings.  When 
this  orderly  migration  occurs,  no  other  anomalous  behavior  ensues,  i.e.,  no  ESF,  no 
TEP,  and  no  scintillation. 

On  other  evenings,  some  mechanism  (possibly  gravity  waves)  disturbs  the 
orderly  migration  of  electrons  back  toward  the  equator.  A  perturbation  can  cause 
turbulence  to  erupt  within  the  anomalies.  This  evolution  can  be  understood  by  first 
considering  the  undisturbed  situation  shown  in  Figure  6.  As  seen,  the  viewpoint  is 
looking  north  at  the  juncture  of  the  ionosphere  and  the  free-space  region  of  the 
stratosphere.  The  magnetic  field  lines  are  pointing  north  into  the  page  and  the 
electron  density  is  N  electrons/m3.  Gravity  attempts  to  force  the  charged  particles 
to  descend  in  the  direction  of  the  acceleration,  g,  but  the  particles  are  held  in  place 
by  the  magnetic  field.  The  manifestation  of  the  opposing  force  of  B  is  a  current  flow 
generated  in  the  direction  indicated  by  J  in  the  figure.   In  an  electrically  neutral 
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Figure  6.  Current  produced  by  gravity/B  interaction. 


ionosphere  (the  number  of  ions  equals  the  number  of  electrons),  the  current  density 
is  given  by  Chen  (1984)  as 


J  =  N(M  +  m) 


g  x  B 


(7) 


where  m  and  M  are  the  electron  and  ion  masses,  respectively.  Current  flowing  as 
indicated  on  the  figure  creates  an  electric  field  pointing  in  the  same  direction.  The 
resulting  E  x  B  force  opposes  the  g  force  and  keeps  the  charged  particles  from 
falling,  resulting  in  a  steady-state  condition. 

This  steady-state  becomes  unstable  in  the  presence  of  an  external  electric 
field.  In  Figure  7  the  equilibrium  situation  of  Figure  6  has  been  altered  by  the 
presence  of  a  transient  electric  field.  With  positive  alternations  of  the  transient  field, 
the  upward  directed  E  x  B  force  is  augmented  by  AE  x  B,  which  overcomes  the 
force  of  gravity,  causing  the  plasma  to  ascend.  With  negative  alternations,  the 
opposite  effect  occurs,  forcing  the  plasma  in  this  region  downward.  The  vortex 
formed  by  plasma  motion  destabilizes  the  ionosphere  within  the  anomalies,  creating 
a  rising  bubble  which  is  devoid  of  plasma.  The  plasma  which  was  contained  in  the 
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AEXB 
Figure  7.    Transient  electric  field  influence  on  E  x  B 
field. 


plasma-deficient  bubble  is  forced  to  the  top  of  the  bubble.  Initial  stages  of  bubble 
formation  was  captured  by  the  Altair  non-coherent  scatter  radar  at  Kwajalein  Island 
and  is  shown  in  Figure  8.  The  deficient  region  extends  from  the  base  of  the 
ionosphere  near  300  km  to  beyond  500  km. 
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Figure  8.    Initial  bubble  formation  over  Kwajalein  (From 

Tsunoda,  1981) 
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Observations  indicate  an  ascending  of  the  deficient  region  similar  to  that 
observed  when  air  bubbles  rise  within  a  heavier  covering  liquid,  demonstrating  a 
Rayleigh-Taylor  instability.  As  the  plasma-deficient  region  of  the  bubble  rises  into 
the  denser  area,  it  spreads  north  and  south  along  magnetic  field  lines  and  grows 
higher  and  elongates  toward  the  magnetic  poles.  Steep  density  gradients  form  at  the 
boundary  between  the  deficient  area  and  the  surrounding  highly  ionized  area, 
forming  an  elongated  duct.  (Kelley,  1989) 

Mature  ducts  extend  to  locations  near  the  anomalies  (i.e.,  ±20°  ).  Although 
their  dimensions  vary,  an  average  configuration  has  been  observed  to  be  ovate  with 
the  east-west  dimension  approaching  50  km  and  the  vertical  dimension  spanning  200 
km.  The  electron  density  within  the  bubble  itself  is  very  low  but  at  the  ceiling  it  is 
high,  with  a  density  differential  of  three  orders  of  magnitude  (Tsunoda,  1981).  The 
structure  developed  by  this  formation  is  an  ionospheric  duct,  which  can  act  as  a 
waveguide,  that  traps  north-south  traveling  electromagnetic  waves  and  allows  them 
to  propagate  over  great  distances  to  the  opposite  hemisphere.  The  structure  of  the 
duct  ceiling  is  not  homogeneous;  it  is  imbedded  with  tube-like  irregularities  of 
varying  size.  Each  irregularity,  depending  on  its  size  and  density,  bends  and  delays, 
by  differing  amounts,  electromagnetic  signals  passing  through  it,  which  can  cause 
satellite  signal  scintillations.  A  pictorial  representation  of  an  ionospheric  duct  is 
shown  in  Figure  9. 

After  formation,  the  ducts  migrate  eastward  in  the  nighttime  equatorial 
ionosphere,  due  to  prevailing  nighttime  eastwardly  winds  at  F-region  heights  (Fejer, 
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Figure  9.  Top  and  side  views  of  an  ionospheric  duct 

et  al.,  1981).  Duct  formation  is  varied;  single  and  multiple  ducts  can  form,  all 
migrating  eastward.  An  overhead  representation  of  these  migrating  ducts  is  shown 
in  Figure  10.  The  imbedded  irregularities  move  with  the  duct,  producing  different 


Figure  10.    Ionospheric  ducts  moving  eastward  in 
nighttime  sky 

amounts  of  signal  path  bending,  attenuation,  delay,  and  Doppler  frequency  shift  for 
an  electromagnetic  wave  passing  through  the  duct.     This  variation  in  signal 
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amplitude,  phase,  and  frequency  is  called  scintillation.  Although  an  overhead  duct 
may  migrate  eastward  and  no  longer  influence  the  area  of  the  earth  where  it  formed, 
other  ducts  which  may  have  formed  to  the  west  progress  overhead  and  continue  to 
affect  local  observations. 

Propagating  electromagnetic  waves,  incident  on  duct  apertures  at  appropriate 
angles  can  become  trapped  within  the  ducts.  Propagation  within  a  duct  is  analogous 
to  waveguide  propagation.  Ducts  follow  magnetic  lines  of  flux  which  have  large  radii 
of  curvature  (Hargreaves,  1992). 
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III.  REMOTE  OBSERVATION  OF  THE  EQUATORIAL  IONOSPHERE 

Knowledge  of  equatorial  ionospheric  conditions  is  essential  for  both 
conducting  scientific  ionospheric  studies  and  for  updating  real-time  employment. 
For  example,  the  presence  of  an  overhead  duct  indicates  whether  or  not  TEP  is 
possible  or  whether  there  might  be  satellite  communications  scintillation.  There  are 
several  methods  by  which  the  presence  of  an  overhead  duct  can  be  ascertained. 
Regardless  of  the  ionospheric  observation  method  used,  results  may  not  apply 
beyond  a  few  hundred  miles  in  longitude.  Propagation  conditions  may  not  be  the 
same  or  even  similar  at  a  location  distant  from  an  observation  site. 

Remote  ionospheric  observation  can  be  accomplished  by  using  equipment 
based  either  on  the  earth  or  on  satellites.  Satellite-based  observation  equipment 
provides  important  information  about  the  topside  portion  of  the  ionosphere  which 
cannot  be  observed  from  earth.  Since  the  goal  of  this  research  was  to  help 
characterize  the  behavior  of  the  lower  ionosphere,  particularly  with  respect  to 
equatorial  spread-F  and  ducts,  only  earth-based  observation  equipment  was  used. 

Active  earth-based  observation  equipment  includes  coherent  and  non-coherent 
scatter  radars,  and  vertical  and  oblique  incidence  sounders.  Passive  devices  are 
Faraday  rotation  and  time-of-flight  measuring  equipment,  and  atomic  recombination 
airglow  monitors.  The  equipment  used  for  ionospheric  monitoring  during  this 
experiment  consisted  of  an  oblique  sounder,  apparatus  to  measure  time-of-flight  of 
transionospheric  signals,  and  an  airglow  monitor. 
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A.         IONOSPHERIC  SOUNDERS 

Ionospheric  sounders,  also  called  ionosondes,  are  monostatic  or  bistatic.  The 
monostatic  ionosonde  is  configured  to  monitor  the  ionosphere  directly  overhead  in 
the  vertical  incidence  configuration,  or  with  an  obliquely-launched  antenna  elevation 
angle,  to  monitor  the  ionosphere  at  a  distant  point  using  ionospheric  backscatter. 
The  bistatic  ionosonde  is  used  exclusively  to  observe  ionospheric  conditions  over  the 
path  between  the  transmitting  and  receiving  antennas  in  the  oblique  incidence  mode. 

Regardless  of  type,  ionosondes  emit  pulsed,  swept-frequency  signals, 
analogous  to  that  of  a  chirp  radar.  The  time-of-flight  is  measured  from  transmission 
to  reception,  from  which  the  ionospheric  layer  heights  can  be  inferred.  Sweeping  the 
frequency  allows  the  electron  density  of  the  layers  to  be  ascertained. 

The  output  data  from  sounders  is  stored  graphically,  depicting  the  layer 
heights  on  an  ionogram.  Newer,  digital  ionosondes  (digisondes)  store  the  data  on 
magnetic  media  which  can  be  evaluated  using  analysis  software  or  displayed  as  an 
ionogram. 

Although  not  used  in  this  experiment,  the  vertical  incidence  sounder  is 
described  first.  Oblique  incidence  sounder  theory  is  derived  from  that  of  the  vertical 
incidence  sounder. 

1.        Vertical  Incidence  Ionosondes 

Creating  ionograms  using  the  vertical  incidence  ionosonde  (VI)  is  the  oldest 
method  for  monitoring  the  overhead  ionosphere,  and  is  still  widely  used  today.  A 
typical  VI  sweeps  over  a  frequency  range  of  0.1  to  30  MHz  in  a  time  period  of  a  few 
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seconds  to  a  few  minutes  (Davies,  1990).  A  pulse  at  each  frequency  is  transmitted, 
after  which  the  sounder  switches  to  the  receive  mode  for  a  certain  length  of  time. 
If  an  echo  is  received  at  that  frequency,  the  round  trip  time  of  flight  is  measured,  and 
a  computed  layer  height  is  recorded.  The  next  frequency  is  then  transmitted  and  the 
process  repeated. 

Echos  are  created  by  the  electromagnetic  pulse  being  apparently  reflected  from 
an  ionospheric  layer.  Although  reflection  is  the  assumed  mechanism  for  returning 
the  energy,  it  is  actually  returned  through  the  process  of  refraction.  When  waves 
propagate  from  a  medium  of  one  refractive  index  into  a  second  region  with  a 
different  refractive  index,  the  wavefront  is  bent  or  refracted.  Ignoring  collisions  and 
the  magnetic  field,  the  refractive  index  of  the  ionosphere  is  given  by  Davies  (1990) 


as 
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where  fN  is  the  plasma  frequency.  The  plasma  frequency  is  a  quantity  dependent 
upon  physical  constants  and  the  plasma  electron  density,  i.e., 
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where  e  and  m  are  the  charge  and  mass  of  an  electron  and  N  is  the  electron  density 
in  electrons/m3.  These  variables  are  all  constants  except  for  N  so  that 
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As  the  electron  density  increases,  the  plasma  frequency  rises,  and  for  a  given 
frequency,  n  decreases. 

The  refractive  index  for  the  ionosphere  is  therefore  a  function  of  frequency. 
The  range  over  which  the  refractive  index  varies  can  be  divided  into  three  frequency 
dependent  regions:  (1)  f  >  fN,  so  that  0  <  n  ^  1,  (2)  f  =  fN,  and  n  =  0,  and  (3)  f  < 
fN,  yielding  n  <  0.  Application  of  SnelPs  law  to  these  three  cases  reveals  that  in 
region  1  the  wavefront  is  bent  toward  the  denser  medium,  in  region  2  the  wavefront 
is  bent  90  degrees  indicating  a  maximum  height  is  attained  (called  the  height  of 
reflection  hj,  and  in  region  3,  the  refractive  index  is  imaginary,  signifying  an 
evanescent  region.  The  conclusions  to  be  drawn  from  these  observations  are  if  (1) 
f  >  fN,  the  propagating  wavefront  will  be  bent  but  not  returned,  (2)  f  =  fN,  the  wave 
is  reflected  from  the  layer,  and  (3)  f  <  fN,  the  wave  will  not  propagate  within  that 
region. 

The  VI  therefore  operates  by  monitoring  echos  as  the  frequency  is  increased. 
The  frequency  at  which  an  echo  no  longer  returns  is  equal  to  the  plasma  frequency 
for  the  layer,  from  which  N  can  be  computed.  The  height  at  which  this  occurs  is  the 
layer  height.  This  operation  is  conducted  for  each  of  the  reflecting  ionospheric 
layers,  i.e.,  the  E,  Fp  and  F2  layers.  The  plasma  frequencies  for  these  three  layers  are 
named  f0E,  f^,  and  f0F2,  respectively.  An  idealized  daytime  vertical  incidence 
ionogram  is  shown  in  Figure  11.  At  night  the  ionogram  consists  of  just  two  traces, 
as  the  F:  disappears  after  sunset.  Above  the  plasma  frequency  for  each  layer,  the 
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transmitted  signal  does  not  return  an  echo,  but  instead  escapes  to  a  higher  layer  or 


3  4  5  6  7 

Frequency  (MHz) 

Figure  11.  Ideal  vertical  incidence  ionogram  (After  Hunsucker,  1991). 


The  presence  of  ESF  in  a  vertical  ionogram  is  the  result  of  overhead  ducts. 
Therefore,  the  VI  can  be  utilized  to  signal  the  presence  of  an  ionospheric  duct 
overhead.  The  manifestation  of  ESF  on  the  vertical  incidence  ionogram  is  that  the 
reflections  from  the  F  layer  are  not  the  clean,  coherent  echos  as  shown  in  the  figure 
above.  Instead,  the  reflecting  layer  becomes  dispersed  in  altitude,  spreading  the  F- 
layer  return.  An  example  of  an  equatorial  spread-F  vertical  incidence  ionogram  is 
shown  in  Figure  12. 

The  first  frame  of  the  figure  begins  in  the  upper  left  corner  where  four 
ionogram  traces  are  seen.  The  bottom  trace  represents  the  E  layer,  while  the  other 
three  depict  the  F  layer.  The  lowest  F  layer  trace  is  the  actual  F-layer  return  and  the 
other  two  are  round  trip  reflections  from  the  earth.  In  the  first  frame,  the  E  and  first 
F  layer  returns  are  clean,  nearly  ideal  traces. 
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Figure  12.  Equatorial  spread-F  progression.  Beginning  in  upper  left  corner,  then 
down,  second  column,  etc.,  the  times  of  each  frame  are:  2215,  2230,  2245,  2300, 
0000,  0100,  0300,  0330,  0400,  0500,  0515,  0530,  0545,  0600,  0645,  and  0700. 
(From  McNamara,  1991) 


Directly  below  the  first  frame,  fifteen  minutes  later,  the  second  frame  shows 
the  first  F-trace  beginning  to  acquire  some  "fuzz"  or  spread  in  height.  The  F-spread 
can  be  seen  in  the  later  frames  to  progressively  deteriorate  until  morning  when  the 
layers  re-form.  From  these  ESF  returns,  overhead  duct  presence  can  be  easily 
ascertained  due  to  the  spread  on  the  F  traces.  It  will  be  shown  that  ESF  is  readily 
recognizable  from  oblique  incidence  ionograms,  as  well. 
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2.         Oblique  Incidence  Ionosondes 

The  oblique  incidence  ionosonde  (OI)  is  situated  with  its  transmitter  and 
receiver  spatially  separated,  usually  over  great  distances.  Depending  upon  the 
number  of  ionospheric  hops  between  transmitter  and  receiver,  several  ionospheric 
reflection  points  can  contribute  to  the  composition  of  the  ionogram.  To  preclude  the 
ambiguities  caused  by  the  multi-hop  path,  the  circuit  planner  usually  limits  the  path 
length  to  that  of  one  hop,  approximately  4000  km  in  the  mid-latitudes.  When 
sounding  across  the  magnetic  equator,  the  single  hop  distance  can  be  extended 
beyond  4000  km  due  to  chordal  reflection  by  the  anomalies  during  the  day,  and 
trapping  within  an  ionospheric  duct  at  night. 

The  operation  of  the  OI  is  similar  to  that  of  the  VI.  However,  the  transmitter 
and  receiver  of  the  VI  are  co-located  and  operate  from  the  same  precise  clock.  The 
OI  operates  from  two  separate  clocks  which  must  be  synchronized  in  order  to 
accurately  measure  time  of  flight  of  the  propagating  signals.  Modern-day  sounders 
utilize  GPS  timing  to  maintain  synchronicity  between  the  clocks. 

The  OI  ionogram  differs  from  that  of  the  VI  in  two  principal  respects.  First, 
the  plasma  frequency  of  the  layers  is  effectively  increased,  so  that  a  higher  frequency 
can  be  reflected  by  an  ionospheric  layer,  and  second,  a  high  or  Pederson  ray 
propagates  obliquely  (but  not  vertically),  modifying  the  appearance  of  the  oblique 
ionogram.  The  reflection  frequency  increase  will  discussed  first. 

The  angle  formed  at  the  ionospheric  entrace  point  by  the  Poynting  vector  of 
the  propagating  wave  and  the  vertical  is  called  the  incidence  angle,  <j)0.     (The 
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incidence  angle  of  a  VI  wave  is  zero.)  Snell's  law  states  r^  sin  <j>0  =  nx  sin  <)),  =  const. 
The  first  region  is  free  space,  i.e.,  r^  =  1,  and  the  second  region  is  at  the  point  of 
reflection.  Snell's  law  indicates  sin  <f>0  =  n,.  sin  <j>r,  where  the  r  subscript  betokens  the 
reflection  point.  The  incidence  angle  at  the  point  of  reflection  is  <|)r  =  90°  producing 
a  refractive  index  at  reflection 

n  =  sin  4>0.  (11) 

By  substituting  this  value  for  n  into  Eq.  8,  the  frequency  that  reflects  from  the  layer 
is  found  as 

/  =  fN  sec  <j>0.  (12) 

Since  sec  <J)0  >  1  for  oblique  incidence,  the  frequency  that  will  reflect  is  always 

greater  than  that  at  vertical  incidence. 

The  effect  of  higher  frequencies  being  reflected  on  the  oblique  ionogram  is 
that  the  layer  traces  of  Figure  1 1  will  be  pushed  to  the  right,  i.e.,  an  individual  trace 
will  bend  toward  the  vertical  at  a  higher  frequency  as  shown  in  Figure  13.  Also 
shown  in  the  figure  are  the  results  of  the  second  modifier  of  the  OI  ionogram,  the 
Pedersen  ray,  which  is  discussed  next. 

As  a  reflection-bound  ray  with  frequency  f  progresses  obliquely  into  an 
ionospheric  layer,  it  eventually  reaches  the  point  in  the  layer  at  which  f  =  fN  sec  c|>0 
where  it  is  reflected.  If  this  point  is  not  at  the  peak  electron  density  of  the  layer,  the 
plasma  frequency  increases  at  altitudes  higher  than  the  reflection  point.  This  factor 
allows  for  the  possibility  of  a  second  ray  satisfying  the  requirements  of  Eq.  12. 
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Figure  13.  Ideal  oblique  incidence  ionogram. 
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In  Eq.  12,  for  a  fixed  frequency  f,  there  are  two  degrees  of  freedom  which  can 
satisfy  the  equality.  The  first  is  found  by  solving  for  the  incidence  angle  when  fN  is 
fixed  at  a  value  less  than  layer  maximum.  The  other  is  established  by  letting  fN 
assume  the  plasma  frequency  value  at  the  peak  of  the  layer.  The  equality  can  still 
be  satisfied  by  decreasing  the  incidence  angle.  This  second  ray  is  called  the  Pedersen 
ray.  Since  this  ray  ascends  higher  into  the  layer  than  the  first,  the  Pedersen  ray  is 
also  called  a  high  ray.  When  fN  for  the  first  or  low  ray  occurs  at  the  peak  of  the  layer, 
there  is  only  one  ray  for  the  path. 

The  Pedersen  ray  produces  its  own  trace  on  the  OI  ionogram.  For  every 
frequency  transmitted  there  are  two  paths  for  the  signal  to  propagate  from 
transmitter  to  receiver,  producing  two  traces  on  the  ionogram  (Nielson  and  Hagn, 
1964).   The  second  trace  is  not  always  detected  however,  since  it  suffers  greater 
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attenuation  over  its  longer  path  length.  When  the  transmitted  frequency  reaches  the 
peak  plasma  frequency  for  the  layer,  the  two  traces  converge  into  one. 
B.         GLOBAL  POSITIONING  SYSTEM 

Sounders  can  be  used  to  determine  peak  electron  density  and  height.  They 
convey  no  information  concerning  the  total  electron  content  of  the  ionosphere, 
however.  By  measuring  the  time  of  flight  of  transionospheric  signals  using  the 
Global  Positioning  System,  the  ionospheric  total  electron  content,  and  how  it  varies, 
can  be  observed. 

The  Global  Positioning  System  (GPS)  consists  of  a  network  of  twenty-four 
active  semisynchronous  satellites  which  broadcast  continuous,  phase-modulated 
signals  on  two  precise  frequencies.  The  period  of  each  satellite  is  approximately 
twelve  hours.  The  signal  passes  through  the  ionosphere  en  route  from  the  satellite 
to  the  earth-borne  receiver.  Precise  measurements  of  the  time  of  flight  from 
transmitter  to  receiver  can  be  made,  from  which  the  total  electron  content  of  the 
ionosphere  can  be  deduced. 

The  two  signals  broadcast  by  each  satellite,  called  LI  and  L2,  are  transmitted 
on  a  pair  of  coherently  derived  frequencies.  They  are  coherent  in  that  they  are 
integer  multiples  of  a  common  base  frequency,  f0,  where 

f0  =  10.23  MHz.  (13) 

The  frequency  of  LI  is  set  at 

ft  =  qx  f0  =  154/0  =  1575.42  MHz.  (14) 

and  for  L2  is 
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f2  =  q2f0  =  120/0  =  1227.60  MHz.  (15) 

(The  constants  qx  =  154  and  o^  =  120  are  multipliers  of  the  base  frequency.)  Since 
the  ionosphere  is  a  dispersive  medium,  coherency  between  the  two  frequencies 
allows  for  precise  measurement  of  how  the  ionosphere  affects  each  frequency. 

The  primary  function  of  a  GPS  receiver  is  to  measure  the  time  required  for  LI 
and  L2  to  travel  from  a  GPS  satellite  to  the  receiver.  This  is  accomplished  by  the  use 
of  very  stable  local  oscillators  to  measure  the  carrier  phase  of  both  LI  and  L2  with 
great  accuracy.  As  these  "transionospheric"  signals  pass  through  the  ionosphere  en 
route  to  earth,  the  effects  caused  by  the  ionosphere  can  be  detected,  measured,  and 
quantified,  particularly,  the  difference  in  phase  between  LI  and  L2. 

As  shown  in  Section  A,  the  ionospheric  index  of  refraction  is  a  function  of 
frequency  and,  for  frequencies  much  greater  than  the  plasma  frequency,  Eq.  8  can 
be  approximated  by  (Hargreaves,  1992) 
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WN 


n  *  1  -  (w  »  o)J,  (16) 

2o)2 

where  ooN  =  27tfN.  The  value  of  fN  was  given  as  a  function  of  N  in  Eq.  9  so  that 

u2N  =  3182N.  (17) 


Substitution  of  coN2  into  Eq.  16  yields  the  refractive  index 

1591N       .,        40.3N 

n  -  1  -  —  =  1  -  — — ,  (18) 

or  fz 

where  go  is  in  radians/ second  and  f  in  hertz. 
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As  an  electromagnetic  wave  travels  a  distance  C  through  a  medium,  the  phase 
change  through  the  medium  is  Aej>  =  -p£  (where  p  is  the  phase  constant  in  radians 
per  meter),  if  p  is  constant  in  the  medium.  If  p  is  not  constant,  the  phase  must  be 
found  from  integration,  i.e., 

A<J>  =  -  f  p  dlt  =  -J  ^p  eft.  (19) 

After  n  is  inserted  from  Eq.  18,  the  phase  change  after  passing  through  an  ionized 
medium  becomes 


A(j>  =  -  ~  f  n  dH  =  -^ f 
a)  .       1591 


1  _   1591N 
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eft 
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C  +  ^--^   I  N  eft. 
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The  first  term  is  simply  the  phase  delay  that  occurs  after  traveling  a  distance  (!  at  the 
speed  of  light,  c.  The  second  term  represents  the  phase  advance  caused  by  the 
ionized  medium,  where  the  refractive  index  is  less  than  unity.  The  integral  term  is 
called  the  total  electron  content  (TEC)  of  the  ionized  medium, 

TEC  =  f  N  eft,  electrons/ m2.  (21) 

When  the  results  are  combined,  the  phase  change  is 

.  .          to  ,       5.3jc1(T6  rrvr         2n  .       8.44xl(T7  ^r  ,00^ 

Aq>  =  -—  Q  +  TEC  =  - —  {  +  TEC,  (22) 

ceo  X  f 

where  X  =  c/f  and  o  =  2itf. 
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At  the  receiver,  the  phase  can  be  measured  but  the  values  of  Q  and  TEC  which 
contributed  to  it  are  unknown-one  equation  and  two  unknowns.  With  a  two-channel 
GPS  receiver,  a  second  equation  is  added  (the  phase  measurement  on  L2  as  well  as 
LI),  while  the  values  of  H  and  TEC  do  not  change.  For  the  phase  measurements  to 
be  combined  in  a  meaningful  way,  they  are  normalized  to  the  base  frequency  by 
dividing  by  their  frequency  multipliers,  qx  and  c^.  The  normalized  phase 
measurements  of  LI  and  L2  are 

,  2n    „       8.44xl(r7  __         27I/o  „       8.44*10'7  _.       .      ,    _ 

9.  =  - — —  H  + TEC  = Q  +  - TEC,      1  =  1,  2.  (23) 
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The  difference  between  these  normalized  phase  measurements  yields 
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Solving  for  TEC, 


TEC  =  4.44X1017  5*  ^^ 


m 


-  44.4  S4>  7TC  units, 


(25) 


where  1  TEC  unit  is  equivalent  to  1  x  1016  electrons/m2. 

The  TEC  is  a  measure  of  the  number  of  free  electrons  in  the  ionosphere 
between  the  observer  and  the  satellite.  In  a  stratified  ionosphere,  as  is  normally 
assumed  at  mid-latitudes,  the  TEC  remains  virtually  constant  throughout  the  day. 
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There  are  gradual  changes  associated  with  the  amount  of  ionizing  energy  from  the 
sun,  so  that  at  noon  the  TEC  is  expected  to  be  at  its  peak  and  diminished  somewhat 
at  other  times.  In  the  tropical  ionosphere,  stratification  is  not  expected  nor  usually 
observed.  As  ducts  move  overhead,  the  number  of  free  electrons  varies  continuously, 
producing  the  subsequent  changing  of  the  TEC  between  observer  and  satellite.  From 
the  algorithm  described  in  Eq.  25,  the  variation  of  TEC  can  readily  be  derived  from 
measurements  over  time  of  the  phase  of  LI  and  L2  at  the  receiver. 

Once  the  variation  of  TEC  is  determined,  its  value  can  be  substituted  into  Eq. 
22  to  find  the  ensuing  fluctuation  of  the  phase  of  the  received  signal.    This 
fluctuation  is  the  phase  scintillation.   Phase  scintillation  can  prevent  the  receiver 
from  being  able  to  tune  and  lock  on  to  the  incoming  signal. 
C.         AIRGLOW  MEASUREMENT 

The  earth's  atmosphere  is  comprised  primarily  by  nitrogen  and  oxygen  with 
lesser  amounts  of  argon,  carbon  dioxide,  helium,  and  hydrogen.  As  the  altitude 
above  the  earth  varies,  the  partial  pressures  contributed  by  each  component  varies. 
At  different  altitudes,  different  components  become  predominant.  Between  the 
altitudes  of  about  200  -  650  km,  the  predominant  component  is  atomic  oxygen. 
(Hargreaves,  1992) 

When  solar  energy  ionizes  atoms  and  molecules,  the  free  electrons  begin  to 
recombine  with  positively  charged  ions.  Two  types  of  recombination  occur: 
molecular  and  atomic,  examples  of  which  are 


34 


NO  *  +  e  ~  -  N  +  O,         (molecular)  and 

(26) 
0+  +  e "  -  O  +  hv      (atomic). 


The  atomic  recombination  emits  a  photon  with  energy  hv  where  h  is  Planck's 
constant  and  v  is  the  photon  radiation  frequency.  The  molecular  recombination 
occurs  at  a  rate  of  about  1000  times  that  of  the  atomic  recombination  (Kelley,  1989). 
The  second  reaction  is  said  to  be  radiative  as  a  photon  is  emitted. 

After  sunset,  solar  ionization  ceases  and  no  new  ions  are  created.  Molecular 
recombination  occurs  rapidly  leaving  only  the  atomic  ions  available  for  further 
recombinations.  Between  200  and  650  km,  in  the  F-region,  the  predominant 
component  is  atomic  oxygen,  constraining  the  majority  of  the  ions  to  be  0+. 
Recombination  with  0+  is  radiative  and  a  photon  is  released  as  the  recombined  atom 
returns  to  ground  state  energy  level.  The  predominant  energy  level  of  these  photons 
is  such  that  they  radiate  a  red  color  at  630  nm  wavelength  (Hargreaves,  1992).  It 
has  been  determined  that  630  nm  emission  originates  predominantly  between  200 
and  400  km  (Weill,  1966).  Radiative  emissions  in  the  ionosphere,  such  as  these,  are 
called  airglow. 

Airglow  is  measured  by  the  photon  emission  rate  in  Rayleighs,  R,  where 

1R  =  106  photons  cm  '2s  \  (27) 

Emission  rates  are  directly  proportional  to  the  electron  density  (Steiger,  1966)  so 

that  measurement  of  the  airglow  recombination  rate  is  an  indirect  method  of 

measuring  electron  density.    Observation  of  airglow  is  accomplished  by  using  a 
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photometer  which  is  an  optical  device  with  a  large  viewing  angle  but  a  small 
bandwidth.  Typical  photometers  have  bandwidths  of  a  few  tens  of  angstroms. 

To  summarize  the  foregoing  paragraphs,  by  measuring  the  630  nm  airglow, 
the  electron  density  of  the  F-region,  particularly  between  200  and  400  km,  can  be 
ascertained. 

Airglow  photometers  can  be  used  near  the  magnetic  equator  to  monitor  the 
low  latitude  ionosphere.  Since  overhead  ducts  have  a  deficiency  of  electrons  and 
ions,  there  are  fewer  recombinations  than  outside  the  duct.  With  fewer 
recombinations,  the  airglow  intensity  is  diminished  looking  up  into  the  duct.  As  the 
ducts  move  eastward  in  the  nighttime  sky,  a  contrast  can  be  seen  between  the  dimly 
glowing  ducts  and  the  separating  ionosphere.  Airglow  deficient  regions  have  been 
observed  in  the  equatorial  region  to  extend  to  about  1200  km  north-south  and  50- 
200  km  east-west  (Weber,  et  al.,  1978). 
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IV.  ANOMALOUS  RADIOWAVE  PROPAGATION  AT  LOW  LATITUDES 

Propagation  of  radiowaves  which  does  not  follow  accepted  behavior  is 
considered  anomalous.  In  low  latitudes,  non-standard  propagation  caused  by  the 
ionosphere  can  be  traced  to  the  equatorial  anomalies  and  the  instabilities  that  arise 
upon  their  dissipation.  Because  of  the  high  electron  density  of  the  anomalies,  the 
sphere-of-influence  of  the  ionosphere  is  no  longer  limited  to  just  the  High  Frequency 
(HF)  band  as  in  mid-latitudes,  but  extends  into  VHF  and  beyond. 

This  low-latitude  anomalous  behavior  can  be  divided  into  two  domains: 
extended  range  propagation  across  the  magnetic  equator  on  terrestrial  paths,  and 
inordinate  scintillation  on  earth-space  links.  The  first  domain  falls  under  the  general 
heading  of  transequatorial  propagation  while  the  second  appertains  to 
transionospheric  propagation. 
A.         TRANSEQUATORIAL  PROPAGATION 

Sky  wave  propagation  entails  a  radiowave  transmitted  by  an  earth-bound 
antenna  in  the  skyward  direction,  with  subsequent  ionospheric  reflections  before 
being  received  at  a  distant  location.  If  the  wave  experiences  only  one  ionospheric 
reflection  en  route,  the  path  is  termed  single  hop.  Due  to  the  spherical  geometry  of 
the  earth  and  the  ionospheric  layers,  standard  single-hop  paths  are  limited  to 
approximately  4000  km.  The  frequencies  which  are  supported  over  these  long-range 
paths  are  restricted  to  HF. 
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Long-range  single-hop  sky  wave  propagation  in  the  north-south  direction 
across  the  magnetic  equator  can  occur  over  distances  much  greater  than  4000  km. 
The  literature  contains  many  cases  of  established  links  in  excess  of  5000  km  at 
frequencies  into  the  VHF  band.  Some  reports  of  low-UHF  propagation  have  also 
been  documented.  Anomalous  propagation  of  this  genre  is  known  as  transequatorial 
propagation  (TEP). 

Transequatorial  propagation  can  be  categorized  into  two  separate  types.  The 
first  is  a  daytime  phenomenon  involving  the  anomalies,  called  afternoon-TEP.  The 
second  occurs  after  sunset  utilizing  ionospheric  ducts,  called  evening-TEP. 

1.        Afternoon-TEP 

The  anomalies  are  created  by  the  fountain  effect  which  repositions  free 
electrons  created  over  the  magnetic  equator.  This  produces  two  regions  of 
anomalously  high  electron  density.  It  was  shown  in  Chapter  III  that  an  oblique 
incidence  sky  wave  will  reflect  from  the  ionosphere  at  altitudes  where  f  =  fN  sec  4>0. 
The  plasma  frequency  is  proportional  to  the  electron  density  so  that  the  inordinately 
high  electron  density  of  the  anomalies  will  reflect  higher  than  normal  sky  wave 
frequencies. 

The  high  electron  density  of  the  anomalies  explains  how  higher  than  normal 
frequencies  can  be  propagated,  but  it  does  not  justify  the  increase  in  single-hop 
propagation  distance.  The  reason  for  increased  propagation  distance  is  visualized 
by  observing  that  the  radiowave  can  reflect  from  one  anomaly  to  the  other  and 
thence  back  to  earth,  as  is  depicted  in  Figure  14. 
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2.         Evening-TEP 

After  sunset  the  migration  of  the  anomaly  electrons  back  toward  the  magnetic 
equator  can  become  unstable.  As  was  discussed  in  Chapter  II,  this  instability  results 
in  the  formation  of  ionospheric  ducts  which  extend  from  one  equatorial  anomaly  to 
the  other.  The  ceilings  of  the  ducts  follow  the  magnetic  field  lines,  raising  the  tops 
of  the  ducts,  over  the  magnetic  equator,  to  several  hundred  km  in  altitude. 

A  duct  can  act  as  a  waveguide,  trapping  a  properly  oriented  north-south 
traveling  radiowave.  Once  trapped,  the  wave  is  subsequently  guided  to  the  opposite 
hemisphere.  Because  of  the  curvature  of  the  duct,  the  reflections  within  the  duct 
occur  in  a  "whispering  gallery"  mode  whereby  only  the  ceiling  of  the  duct  is  used  for 
reflection.  This  results  in  low  reflection  loss  and  allows  high  signal  strengths  to  be 
received  in  the  opposite  hemisphere.  The  losses  which  do  occur  result  primarily 
from  aspect-dependent  coupling  at  the  apertures  of  the  duct,  discussed  next. 
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For  TEP  there  is  no  direct  free  space  path  between  the  transmitter  and  the 
receiver.  The  losses  which  contribute  to  the  path  loss  of  TEP  are: 

1)  free  space  loss  exterior  to  duct,  LpS; 

2)  coupling  losses  entering  and  exiting  the  duct,  1^;  and 

3)  reflection  losses  within  the  duct,  Lr. 

The  total  loss  is  Lj.(dB)  =  LpS(dB)  +  Lc(dB)  +  L^CdB).  These  loss  locations  are 
illustrated  in  Figure  15. 

Duct 


Figure  15.  Path  losses  from  transmitter  to  receiver  via  ionospheric  duct. 

The  free  space  spreading  loss  represented  by  locations  1  are  easily  computed 
if  the  heights  of  the  duct  apertures  are  known;  the  loss  at  positions  2  is  a  function  of 
the  difference  between  the  direction'  of  the  Poynting  vector  of  the  propagating  wave 
and  the  axis  of  the  duct;  and  the  losses  shown  at  position  3  are  normally  quite  low 
since  the  ceiling  electron  density  is  high.  Loss  factors  1  and  3  are  fairly  static  with 
low  variation,  leaving  the  variations  in  the  path  loss  mainly  a  result  of  the  coupling 
losses  of  loss  factors  2. 
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The  coupling  losses  which  result  as  a  difference  in  the  wavefront  propagation 
directions  and  the  axes  of  the  duct  of  factors  2  are  dependent  upon  several  different 
variables.  Among  these  are  the  height  of  the  duct  aperture,  the  height  of  the  apex 
of  the  duct  (which  follows  the  magnetic  flux  lines),  and  the  relative  positions  of  the 
transmitter  and  receiver  sites  on  the  ground  with  respect  to  their  optimum  ground 
positions,  addressed  in  the  following  paragraph.. 

These  physical  height  parameters  of  the  duct  are  not  known  a  priori  but  the 
ranges  of  movement  are  generally  known  and  can  be  statistically  modeled.  The 
optimum  positions  on  earth  for  signal  transmission  and  reception  can  be  determined 
by  using  the  mean  values  of  the  heights.  Studies  of  these  positions  have  been  made 
both  experimentally  and  theoretically  (Heron,  1981;  Piatt  and  Dyson,  1989),  and 
the  optimum  north-south  positions  have  been  found  to  be  approximately  ±18° 
magnetic  latitude.  The  optimum  east-west  positions  are  for  the  magnetic  longitudes 
of  the  two  sites  to  the  same. 

Piatt  and  Dyson,  in  their  theoretical  waveguide  modeling,  show  a  span  of 
reception  locations  on  the  earth  varying  from  14°  to  26°  magnetic  latitude  and  ±2° 
from  the  transmitter  magnetic  longitude.  The  further  the  reception  position  is  from 
the  optimum  position,  the  lower  the  magnitude  of  the  received  signal  will  be  due  to 
coupling  losses. 

When  scattered  worldwide  reports  of  the  transmit  and  receive  locations  are 
considered,  the  latitudinal  limits  of  the  Piatt  and  Dyson  model  are  seen  to  be 
reasonable.   However,  Nielson  (1969)  successfully  established  a  TEP  path  with  a 
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difference  in  longitude  of  6°.    This  indicates  that  the  theoretical  limit  of  ±2°  in 
longitude  difference  is  too  conservative,  and  that  the  limit  should  be  at  least  ±6°. 
B.         TRANSIONOSPHERIC  PROPAGATION 

The  anomalous  propagation  of  TEP  discussed  above  has  beneficial 
ramifications  in  that  the  propagation  distance  and  frequency  range  have  been 
enhanced.  The  transionospheric  effects  are  not  beneficial  and  can  be  disastrous.  The 
effect  upon  radiowaves  which  pass  through  an  ionospheric  duct  is  scintillation  in 
phase,  amplitude,  and  frequency,  sometimes  severe.  Intense  scintillation  has  been 
observed  on  radiowaves  passing  through  the  low  latitude  ionosphere  at  frequencies 
as  high  as  7  GHz  (Hargreaves,  1992). 

Scintillation  causes  problems  for  all  radio  receivers.  It  prevents  receiver  lock- 
on,  demodulation,  and  error-free  reception,  etc.  In  this  report,  scintillation  effects 
upon  geolocation  will  be  emphasized  and  demonstrated.  Geolocation  methods 
include  Angle  of  Arrival  (AOA),  Time  Difference  of  Arrival  (TDOA),  and  Frequency 
Difference  of  Arrival  (FDOA). 

The  cause  of  the  low-latitude  scintillation  was  discussed  in  Chapter  II,  where 
it  was  shown  that  the  duct  ceilings  are  imbedded  with  myriad-sized  tube-like 
enhancements  of  electron  density  irregularities.  Electromagnetic  waves  passing 
through  these  irregularities  are  bent  and  delayed. 

Quantifying  the  effects  upon  transionospheric  radiowaves  can  be  achieved  by 
utilizing  the  total  electron  content  explained  in  Chapter  III.  The  TEC  is  found  using 
GPS  receivers,  and  the  following  scintillation-induced  quantities  can  be  derived. 
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1.  Time  Delay 

The  differential  phase  change  due  to  propagation  through  an  ionized  region 
was  found,  in  Eq.  22  in  Chapter  III,  to  be 

A  ,  8.44xl(T7  rrvr 

A4>™  =  " **£  (28) 

Changes  in  phase  and  time  can  be  related  as 

A<$>  =  2nfAt,  (29) 

and  the  time  delay  associated  with  the  phase  change  of  Eq.  28  is 

A,       8.44x10  7  _       1.34xl0"7  _- 

At  =  TEC  =  TEC.  (3rn 

(2*0  f  f>  C     } 

This  time  delay  has  the  potential  to  degrade  a  TDOA  geolocation  solution. 

2.  Doppler  Shift 

A  frequency  change  can  be  induced  into  a  signal  passing  through  an  ionized 
region.  The  change  in  instantaneous  radian  frequency,  as  a  function  of  time,  is  found 
by  differentiating  the  signal  phase  with  respect  to  time,  i.e.,  Ago  =  A<}>/At.  The  phase 
of  the  transionospheric  signal  was  given  in  Eq.  22,  and  by  differentiating, 


dt{         f 


(31) 


The  change  in  cyclic  frequency  is  found  as 


=  8.44X10-7    d_  _   1.34E10-7  A  iTEQ  Hz.  (32) 

2rt/        dt  f  dt 


Frequency  changes  have  deleterious  effects  on  FDOA  systems. 
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3.         Angular  Deviation 

By  definition,  a  wavefront  undergoing  refraction  changes  propagation 
direction.  Since  the  refracting  ducts  are  in  motion  and  comprised  of  random  shapes 
and  electron  densities,  the  modifications  in  direction  of  the  waves  passing  through 
them  are  difficult  to  ascertain.  If  the  wave  approaches  the  tube  from  directly  above, 
the  angular  change  in  direction  can  be  estimated  from  wedge  refraction  theory 
(Davies,  1990)  as 

c  8.44xl0~7    d    rrm^       40.3    d    rm,^ 

a  = -^- —  (TEC)  =  -^-  —  (TEC),  (33) 

2n       f2  dx  f2     dx  K     J 

where  d/dx  (TEC)  is  the  change  in  the  TEC  in  the  horizontal  direction.  Angular 
direction  changes  can  destroy  the  geolocation  solution  for  AOA  systems. 
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V.  THE  EXPERIMENTAL  DESIGN 

There  are  two  distinct,  independent  manifestations  of  the  presence  of  an 
ionospheric  duct:  (1)  the  trapping  of  a  propagating  radiowave  within  the  duct  and 
subsequent  shunting  to  the  opposite  hemisphere  and  (2)  the  scintillations  induced 
upon  waves  passing  through  the  duct.    These  are  shown  in  Figure  16.    The  two 
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Figure  16.     Terrestrial  signal  trapping  and  space  signal 
scintillation. 


events  are  temporally  and  geographically  correlated  and  the  presence  of  one  implies 
a  probability  of  the  other. 

To  exploit  this  correlation,  first  a  terrestrial  radio  link  was  established 
between  the  Hawaiian  Islands  in  the  north  Pacific  and  the  Cook  Islands  in  the  south 
Pacific.  The  geography  is  shown  in  Figure  17. 
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Figure  17.  Central  Pacific  showing  Hawaii,  Cook  Islands, 
and  Christmas  Island. 


The  transmitted  signal  was  commercial  television  station  KHON  Channel  2  in 
Honolulu,  between  54  and  60  MHz.  A  receiver  was  placed  in  Rarotonga  in  the  Cook 
Islands  with  a  recorder  calibrated  to  record  received  signal  power  of  the  video  carrier 
at  55  MHz. 

Concurrently,  a  Global  Positioning  System  receiver  was  collocated  with  the 
receiver  in  Rarotonga  to  monitor  the  total  electron  content  of  the  ionosphere  of  the 
southern  Pacific.  The  Jet  Propulsion  Lab  (JPL)  maintains  a  GPS-TEC  receiver  in  the 
Hawaiian  Islands.  The  data  from  this  receiver  (and  many  others  as  well)  is  posted 
to  their  web  page.  This  data  was  utilized  as  an  ionospheric  observation  of  the 
northern  Pacific.  Since  the  viewing  region  of  these  two  GPS  receivers  did  not 
overlap,  a  third  GPS  receiver  was  placed  near  the  equator  on  Christmas  Island.  The 
viewing  region  of  this  third  receiver  overlapped  the  regions  of  the  first  two. 

46 


The  concurrent  monitoring  of  both  a  TEP  signal  via  radios,  and  of  the 
ionosphere  using  GPS,  provides  information  about  the  nighttime  low-latitude 
ionosphere  and  its  effects  on  radiowave  propagation.  The  addition  of  two 
complementary  pieces  of  monitoring  equipment  for  a  short  time  period  enhanced  ±is 
monitoring  capability.  The  first  was  collocating  an  airglow  monitor  with  the  GPS 
receiver  to  observe  the  nighttime  equatorial  ionosphere  above  Christmas  Island.  The 
other  was  the  installation  of  an  oblique  ionosonde  over  the  same  path  as  the  TEP 
circuit  but  directed  in  the  opposite  direction.  The  sounder  not  only  supplies  data 
concerning  the  bi-directionality  of  the  path,  but  delivers  daytime  propagation  path 
characteristics  as  well. 

The  purpose  of  this  experiment  was  threefold:  1)  determine  factors  which 
control  TEP,  either  deterministically  or  statistically,  2)  determine  if  scintillation  can 
be  correlated  with  the  occurrence  of  TEP,  and  3)  develop  predictive  capability  of  TEP 
and  temporal  and  spatial  occurrence  of  scintillation. 

The  remainder  of  this  chapter  outlines  the  criteria  used  in  the  selection  of  the 
sites  for  the  different  equipment  locations. 
A.         TRANSMITTER  AND  RECEIVER  SITES 

The  equatorial  ionosphere  should  first  be  modeled  based  on  an  ideal,  benign 
environment.  Only  by  modeling  thus  can  any  predictability  of  the  occurrence  of  ESF 
or  ionospheric  ducts  can  be  expected.  There  is  no  location  on  earth  which  provides 
a  completely  benign  environment,  but  candidate  areas  should  be  free  of  radical 
temperature,  wind,  and  magnetic  deviation  changes. 
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In  addition  to  a  moderate  environment,  the  transmit  and  receive  locations 
which  are  to  be  used  in  that  environment  must  be  within  latitudinally  limited  regions 
where  known  TEP  originates  and  terminates,  i.e,  approximately  15°  -  25°  magnetic 
latitude. 

A  third  useful  feature  is  the  ability  to  monitor  the  local  ionosphere  at  the 
transmitter  and  receiver  sites,  and  at  mid-point  of  the  path,  the  magnetic  equator. 
Equatorial  monitoring  should  be  at  or  near  the  sub-propagation  point  of  the  TEP 
signal. 

The  sites  selected  meet  all  the  above  requirements.  The  transmitter  was 
located  on  Oahu  in  the  Hawaiian  Islands  and  the  receiver  was  on  Rarotonga  in  the 
Cook  Islands.  The  first  requirement  is  met  since  the  propagation  path  is  entirely  over 
the  Pacific  Ocean  with  no  mountains  to  alter  the  neutral  winds,  the  magnetic  equator 
is  very  near  the  geographic  equator,  and  the  magnetic  declination  is  about  11°.  The 
second  requirement  is  met  since  TEP  observations  have  been  made  on  this  path  for 
30  years.  The  third  requirement  is  met  by  Christmas  Island  situated  at  2°  north  of 
the  equator,  very  near  the  sub-propagation  path. 
B.         EQUIPMENT  EMPLACEMENT 

1.         Oahu 

The  selected  VHF  transmitter,  TV  Channel  2,  KHON  has  an  omnidirectional 
antenna  with  an  ERP  of  100  kW  located  at  21.32N,  157.85W  geographic,  21.45N, 
91.1 1W  magnetic.    The  station  operates  23  hours  per  day,  with  one  hour  off  for 
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maintenance.  The  maintenance  hour  is  always  after  midnight  and  therefore  does  not 
pose  insurmountable  difficulties  in  observing  nighttime  TEP. 

During  September  1995  the  receive  portion  of  an  oblique  sounder  was 
emplaced  on  Oahu.  The  site  location  was  at  the  Navy's  NCTAMS  on  Wahiawa.  The 
sounder  itself  is  the  property  of  and  on  loan  from  the  Defence  Science  and 
Technology  Organisation  (DSTO)  of  Adelaide,  Australia  (analogous  to  an  arm  of 
DOD).  The  frequency  of  the  sounder  transmitter  swept  from  2  to  70  MHz  which  was 
received  on  three  different  antennas.  The  first  receiving  antenna  for  the  frequency 
range  2-32  MHz  was  a  rosette,  already  in  place  at  Wahiawa,  which  was  not  being 
used  for  other  purposes.  The  second,  which  was  used  for  32  -  55  MHz,  was  a  log 
periodic  array  (LP A),  property  of  DSTO.  The  last  antenna,  an  off-the-shelf  Winegard 
TV  yagi,  was  used  for  55  -  70  MHz.  The  clock  was  maintained  precisely  using  GPS. 
The  received  frequency  and  time  of  flight  were  stored  as  digital  ionograms  on  the 
receiving  suite's  computer  hard  drive. 

2.        Rarotonga 

Two  VHF  receivers  were  located  at  the  receive  site  in  Rarotonga,  at  the 
laboratory  of  the  Cook  Islands  Scientific  Officer,  Mr.  Stuart  Kingan.  The  lab  was 
situated  at  21.22S,  159.74W  geographic,  20.53S,  84.64W  magnetic,  for  a 
transmitter-receiver  separation  distance  of  approximately  4700  km.  The  magnetic 
longitude  separation  between  transmitter  and  receiver  was  6.47  degrees.  The 
receivers  were  tuned  to  the  carrier  frequency  of  the  TV  transmitter  from  KHON 
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Channel  2.  Each  receiver  was  fitted  with  a  low-band  TV  Winegard  LP-Yagi  antenna 
with  pre-amps,  one  polarized  horizontally  and  the  other  vertically.  Calibrated  AGC 
voltage  of  each  receiver  was  fed  to  a  dual  trace  chart  recorder.  The  signal  strength 
was  derived  by  conversion  of  the  calibrated  AGC  whenever  the  ducts  exist  and  allow 
TEP  from  the  Oahu  transmitter. 

During  September  1995,  the  transmitter  for  the  oblique  sounder  was 
collocated  with  the  VHF  receivers  in  Rarotonga.  As  with  the  receiver,  the  clock  was 
maintained  with  GPS  timing.  The  swept  sounder  frequency  signal  was  transmitted 
from  two  antennas,  twenty-four  hours  per  day.  The  first  antenna,  used  for  2  -  30 
MHz,  was  an  inverted  Vee.  The  second,  used  from  30  to  70  MHz,  was  a  terminated 
rhombic. 

A  GPS  TEC  system  was  collocated  at  the  receiver  site.  The  GPS  equipment 
was  a  dual  receiver,  12  channel,  Z-12  Ashtech.  With  two  receivers,  both  LI  and  L2 
can  be  received  simultaneously  for  up  to  12  different  satellite  transmitters.  The  TEC 
measurements  were  made  and  recorded  continuously. 

3.        Christmas  Island 

A  GPS  TEC  Z-12  Ashtech  system  was  emplaced  at  the  NOAA  site,  2.01N, 
157.4W  geographic,  2.64N,  86.86W  magnetic,  on  Christmas  Island.  The  TEC  was 
monitored  and  recorded  continuously.  For  a  short  period  in  September  1995,  an 
airglow  monitoring  system  was  erected  and  overhead  measurements  made  and 
recorded  nighdy  by  Dr.  Michael  Taylor  from  Utah  State  University. 
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C.         PREVIOUS  DATA 

An  added  bonus  to  this  particular  path  is  that  Mr.  Kingan  has  recorded  the 
start  and  stop  times  of  the  reception  of  this  same  TV  signal  at  Rarotonga  in  the  past. 
His  data  are  continuous  from  April  1966  through  December  1976.  He  has  graciously 
provided  this  data  which  has  been  incorporated  into  the  analysis. 
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VI.  THE  TRANSEQUATOR1AL  VHF  EXPERIMENT 

The  data  collected  during  the  experiment  outlined  in  Chapter  V  will  be 
analyzed,  and  when  appropriate,  prediction  algorithms  developed.  The  data 
analyzed  in  this  chapter  are  the  VHF  reception  measurements,  which  were  collected 
at  Raro tonga  from  April  1966  to  December  1976  and  from  March  1995  to  March 
1996.  The  early  data  consists  of  the  start  and  stop  times  of  the  reception  of  TV 
Channel  2,  at  55  MHz,  from  Honolulu.  The  latter  data  are  a  continuous  monitoring 
of  the  signal  strength  received  from  Channel  2. 

In  Chapter  V  it  was  noted  that  two  receivers  were  used  to  monitor  Channel 
2,  one  with  vertical  and  the  other  with  horizontal  polarization.  As  expected,  due  to 
the  polarization  becoming  random  during  ionospheric  propagation,  neither 
polarization  delivered  a  power  level  statistically  greater  than  the  other.  In  October 
1995,  one  of  the  receivers  was  re-tuned  to  monitor  other  frequencies  in  addition  to 
Channel  2.  Successful  reception  of  frequencies  up  to  and  including  Channel  11,  at 
200  MHz,  from  Honolulu  has  been  observed  even  though  the  receiving  antenna's 
bandwidth  extends  only  to  Channel  6,  at  about  90  MHz. 

Analysis  of  the  data  is  concentrated  on  those  aspects  of  the  communications 
medium  necessary  for  the  existence  of  a  circuit.  Predictions  of  the  probabilities  of 
a  circuit  existing,  the  starting  and  stopping  times,  the  frequencies  which  will 
propagate,  and  the  path  loss  are  the  goals.  Several  of  these  factors  are  directly 
controlled  by  the  relative  position  of  the  sun  and  by  the  strength  of  the  solar  flux. 
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A.         CIRCUIT  ONSET  TIME 

The  circuit  commencement  times  for  the  continuous  period  from  April  1966 
to  December  1976  are  plotted  in  Figure  18.     Path  onset  times  are  seen  to  exhibit 
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Figure  18.  Circuit  commencement  times. 

some  form  and  regularity  but  are  also  characterized  by  randomness,  as  are  all 
ionospherically  influenced  signals.  The  zero  readings  at  the  bottom  of  the  graph 
show  days  on  which  the  path  did  not  materialize.  The  first  step  in  the  analysis  of  the 
circuit  commencement  signal  is  ascertaining  correlation  between  path  onset  time  and 
a  physical  parameter,  i.e.,  determining  if  the  path  onset  time  is  the  dependent 
variable  based  on  a  physical  parameter  which  is  an  independent  value.  One 
candidate  independent  variable  is  daily  solar  flux  values  for  each  day  in  the  figure. 
Those  solar  flux  readings  are  shown  in  Figure  19. 
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Figure  19.  Solar  flux. 
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An  effective  method  of  determining  correlation  between  path  commencement 
times  of  Figure  18  and  the  solar  flux  of  Figure  19  is  to  compute  the  cross  correlation 
between  the  two.  The  cross  correlation  between  the  path  onset  time  and  the  solar 
flux  is  shown  Figure  20.  Notice  that  there  is  a  good  correlation  of  about  0.65.  This 
correlation  is  strong  enough  to  establish  solar  flux  as  a  reliable  partial  predictor  of 
path  onset  time. 

A  second  candidate  as  an  independent  variable  controlling  onset  time  of  the 
circuit  is  the  solar  zenith  angle,  or  more  specifically  the  cosine  of  the  solar  zenith 
angle.  The  cosine  of  the  solar  zenith  angle  for  the  same  days  as  above  is  plotted  in 
Figure  21. 
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Figure  20.  Normalized  cross  correlation  of  onset  time  and  solar  flux. 
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Figure  21.  Cosine  solar  zenith  angle. 
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The  cross  correlation  function  of  the  path  onset  time  with  the  cosine  of  the 
solar  zenith  angle  was  computed  and  is  shown  in  Figure  22.  Again  a  correlation  of 
about  0.65  is  seen  between  the  two  signals.  This  also  indicates  a  good  candidate  as 
an  independent  variable  with  partial  control  of  circuit  commencement  time. 
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Figure  22.  Normalized  cross  correlation  of  onset  time  and  cosine  solar 
zenith  angle. 


If  days  which  have  no  propagation  are  removed  from  Figure  18,  i.e.,  remove 
the  zeros,  the  cross  correlations  between  circuit  onset  times,  given  that  the  circuit 
exists,  and  the  two  independent  variables  can  be  computed.  These  are  shown  in 
Figures  23  and  24. 
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Figure  23.  Normalized  cross  correlation  of  onset  time  and  solar  flux, 
given  that  circuit  exists. 
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Figure  24.   Normalized  crossed  correlation  of  onset  time  and  cosine 
solar  zenith  angle,  given  that  circuit  exists. 
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In  these  two  figures  cross  correlations  have  increased  to  above  0.95  indicating 
very  strong  correlations  between  circuit  commencement  time  and  each  of  the  other 
two  parameters.  From  this  two  things  are  learned: 

(1)  These  two  candidates  are  indeed  responsible  for  at  least  part  of  the 
response  of  circuit  onset  time,  and 

(2)  The  zeros,  i.e.,  days  when  there  is  no  propagation,  appear  to  be 
controlled  by  one  or  more  other  factors  since  the  correlation  factor  increased  with 
the  removal  of  non-propagating  days. 

The  signals  on  those  days  where  propagation  does  occur  were  smoothed  with 
low-pass  filtering.  Non-causal  finite-impulse-response  filtering  of  50  delay  elements 
was  used  to  obtain  filtered  circuit  onset  time  and  solar  flux.  The  cosine  solar  zenith 
angle  was  also  filtered  to  remove  the  discontinuities  resulting  from  removal  of  the 
non-propagating  days.  These  filtered  factors  are  shown  in  Figures  25,  26,  and  27. 
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Figure  25.  Filtered  onset  times. 
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Figure  26.  Filtered  solar  flux. 
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Figure  27.  Filtered  cosine  solar  zenith  angle. 
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The  filtered  data,  which  emphasize  the  desired  qualities  and  de-emphasize 
those  which  are  undesired,  are  plotted  against  each  other  to  see  how  they  relate  to 
one  another.  The  dependence  of  circuit  commencement  time  upon  solar  flux  is 
obtained  by  plotting  the  filtered  onset  time  against  the  filtered  solar  flux.  The  result 
is  shown  in  Figure  28. 
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Figure  28.  Filtered  onset  time  versus  filtered  solar  flux. 


The  relationship  between  the  filtered  onset  time  and  the  filtered  flux  appears  to  have 
an  inverse  relationship  from  Figure  28.  This  is  also  apparent  by  comparing  the 
general  waveform  shape  of  Figures  25  and  26.  A  good  fit,  in  the  least  squares  sense, 
is  obtained  by  the  quadratic  relationship 

OT(f)  =  0.0000181  (fix)2  -  0.0123(flx)  +  9.46      (GMT),  (34) 

where  OT(f)  is  the  onset  time  response  to  fix,  which  is  the  10.7  cm  solar  flux. 
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Similarly  the  filtered  onset  time  can  be  plotted  against  the  filtered  cosine  solar 
zenith  angle.  In  Figure  29  the  relationship  between  the  two  appears  almost  linear. 
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Figure  29.  Filtered  onset  time  versus  filtered  cosine  solar  zenith  angle. 


A  good  least  squares  fit  is  obtained  by  the  linear  relationship 


OT(x)  =  18.61  -  10.76(cos(x))      (GMT). 


(35) 


The  individual  contributions  of  both  the  zenith  angle  and  the  solar  flux,  and 
the  interaction  between  them  can  be  obtained  by  multiplying  together  Equations  34 
and  35,  i.e., 


OT(fx)  =  OT(f)  *  OT(x). 


(36) 
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The  results  of  the  multiplication  can  then  be  least  squares  fitted  to  the  filtered  onset 
time  to  obtain 


OT  =  0.011  OT(fx)2       1.334  OT(fx)   +  48.06      (GMT). 


(37) 


If  this  equation  is  a  good  fit  to  the  observed  onset  times,  plotting  it  against  the 
data  should  yield  a  straight  line  at  45°.  In  Figure  30  the  filtered  onset  time  is  plotted 
on  the  y  axis  and  OT  from  Equation  37  is  plotted  on  the  x  axis.  Clearly,  the 
prediction  equation  is  valid  since  the  plotted  result  approaches  a  45°  relationship  in 
the  least  squares  sense. 
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Figure  30.  Actual  onset  time  versus  predicted  onset  time. 


Filtered    and    predicted    onset   times   are   plotted    together   against   the 
propagating  days  in  Figure  31.  There  are  a  few  deviations  from  the  observations  in 
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the  predictions,  but  overall  the  data  fit  very  well.  The  observed  data  are  shown  in 
blue  and  the  predicted  values  in  red. 
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Figure  31.  Observed  and  predicted  onset  times. 

Although  the  analysis  of  circuit  onset  time  has  been  conducted  using  GMT 
rather  than  local  time  for  convenience,  the  onset  time  is  actually  a  locally-controlled 
event.  When  predicted  onset  times  are  converted  from  GMT  to  local  time, 


Local  Time  =  GMT  ♦   jgggg* 

15 


(38) 


For  the  equatorial-path  crossing-longitude  of  this  data  set  of  201  °,  the  resulting  local 
time  is  GMT  +  13.4.  With  this  adjustment  Eq.  37  can  be  modified  as 


OT  -  0.011  OT(fxf       1.334  OT(f%)  +  61.46      (local). 


(39) 
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B.         CIRCUIT  CESSATION  TIME 

The  circuit  filtered  cessation  times  for  the  same  continuous  period  from  April 
1966  to  December  1976  are  plotted  in  Figure  32. 
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Figure  32.  Filtered  circuit  cessation  times. 


The  techniques  used  in  the  last  section  to  predict  circuit  onset  times  are  used 
to  determine  a  prediction  model  for  the  circuit  cessation  times  of  Figure  32. 
Independent  variables  are  found  to  be  the  10.7  cm  solar  flux  and  the  sine  of  the  solar 
zenith  angle.  An  exponential  relationship  is  found  which  approximates  the  circuit 
cessation  time  as 


/Zx*sin(x)+10 

CT  =  e        10°         +  10      {GMT). 


(40) 
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A  plot  of  the  predicted  value  of  the  cessation  time  against  the  observed  cessation 
time  is  expected  to  show  a  45°  angle.  This  plot  is  shown  in  Figure  33. 
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Figure  33.  Observed  versus  predicted  circuit  cessation  time. 


In  the  least  squares  sense,  there  is  a  45°  relationship.  A  plot  of  the  predicted  and 
observed  cessation  times  versus  the  days  of  propagation  is  shown  in  Figure  34,  with 
the  observed  cessation  time  shown  in  blue  and  the  predicted  shown  in  red.  There 
is  a  reasonable  match  between  predicted  and  observed  cessation  times. 
By  using  Eq.  38,  Eq.  40  can  be  modified  to  localize  it  as 


CT  =  e 


/k*sin(x)+10 
100 


+  23.4      (local). 


(41) 
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Figure  34.  Comparison  of  observed  and  predicted  cessation  times. 


C.         CIRCUIT  PROBABILITY 

With  the  ability  to  predict  the  starting  and  ending  times  of  the  propagation 
circuit  for  those  days  where  the  circuit  exists  now  available,  a  measure  of  the 
probability  that  the  path  will  actually  exist  is  needed.  The  independent  variables  are 
again  the  relative  position  of  the  sun  and  the  amount  of  ionizing  energy  available  to 
create  electrons  for  the  anomalies  and  the  ensuing  propagation  duct.  In  the  previous 
two  analyses  the  solar  zenith  angle  was  seen  to  be  a  controlling  variable.  In  this 
analysis  of  the  frequency  of  circuit  occurrence,  the  probability  function  will  depend 
more  heavily  on  the  position  of  the  setting  sun  relative  to  the  magnetic  equator.  This 
will  be  defined  in  this  dissertation  as  the  solar  magnetic  angle  or  SMA.  A  method  to 
compute  the  SMA  is  given  in  Appendix  C. 
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The  specific  times  of  occurrence,  for  example  as  shown  in  Figure  18,  are  not 
required  to  determine  if  a  circuit  exists  or  not.  Instead  all  the  starting  times  of  Figure 
18  can  be  converted  to  ones,  i.e.,  the  existence  of  a  path  equates  to  a  one  and  the 
nonexistence  of  a  path  remains  at  zero.  Groups  of  sequential  ones  will  then  indicate 
a  higher  probability  of  occurrence  than  a  single  one  or  a  zero.  If  the  data  are  non- 
causal  filtered  to  provide  moving  averages  of  the  number  of  sequential  ones,  the 
result  is  a  moving  probability  function.  To  obtain  a  finer  scale  than  when  looking  for 
trends  as  in  the  previous  analyses,  only  twenty  delay  elements  are  used.  Figure  35 
shows  the  result  of  converting  the  data  of  Figure  18  to  ones  and  zeros  and  the  non- 
causal  filtering. 
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Figure  35.  Twenty-one  day  moving  average  TEP  circuit  probability. 
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With  knowledge  that  the  circuit  has  the  highest  probability  of  occurrence 
when  the  SMA  is  near  zero,  i.e.,  when  the  sun  sets  on  the  magnetic  equator,  an 
overlaid  normalized  plot  of  the  cosine  of  the  SMA,  scaled  to  the  range  0-1,  and  the 
21 -day  moving  average  circuit  probability,  is  produced  as  shown  in  Figure  36. 
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Figure  36.  Comparison  of  21 -day  moving  average  circuit  probability 
and  cos(SMA).  Indexes  on  plot  indicate  SMA  positions:  1)  Winter 
solstice,  2)  Solar  south-to-north  magnetic  equator  crossing,  3)  Summer 
Solstice,  and  4)  Solar  north-to-south  magnetic  equator  crossing. 


The  cos(SMA)  (in  red)  is  in  phase  with  the  fluctuations  of  the  amplitude  of 
the  probability  plot  (in  blue).  These  magnitudes  are  not  identical,  indicating  that 
there  are  other  controlling  variables.  The  overall  shape  of  the  probability  plot  is 
similar  to  that  of  the  solar  flux  over  the  same  period,  from  Figure  19.  A  good  fit  to 
the  envelope  is  formed  as  a  function  of  the  solar  flux  as 
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fix 
Envelope  =  — —  +  0.6. 
480 


(42) 


If  the  envelope  is  multiplied  by  the  cosine  of  the  SMA  shown  in  red  in  Figure  36,  a 
better  fit  to  the  circuit  probability  shown  in  blue  is  obtained.  The  result  of  this 
multiplication  is  shown  in  Figure  37  with  the  frequency  of  occurrence  signal  again 
shown  in  blue  and  the  modeled  signal  shown  in  red. 
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Figure  37.  21-day  moving  average  circuit  probability  and  cos(SMA) 
multiplied  by  solar  flux  envelope. 


Although  this  model  approaches  the  actual  probability,  another  controlling 
variable  must  be  accounted  for.  Each  cycle  of  the  cosine  of  the  SMA  has  two  peaks, 
and  the  probability  plot  similarly  has  two  peaks  which  are  in  phase  with  the 
cos(SMA)  peaks.  The  first  peak  of  the  probability  shown  in  blue  is  much  lower  than 
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the  second  peak.  This  may  be  due  to  the  delay  in  "charging"  the  local  nighttime 
ionosphere.  Similar  to  a  capacitor  where  the  voltage  lags  the  applied  current,  the 
number  of  free  electrons  in  the  nighttime  ionosphere  is  dependent  upon  the  number 
of  free  electrons  during  the  last  night  and  the  night  before.  Even  though  the  sun  is 
directly  overhead  during  the  day,  there  are  not  as  many  free  electrons  when  the  sun 
approaches  from  the  opposite  hemisphere  as  when  it  approaches  from  the  same 
hemisphere.  Since  the  data  set  has  an  1 1  °  magnetic  declination  angle,  the  sunset 
side  of  the  magnetic  equator  makes  an  angle  of  1 1  °  north  of  the  geographic  equator. 
When  the  sun  passes  the  magnetic  equator  heading  north  [the  first  peak  in  the 
cosine (SM A)],  there  are  relatively  few  free  electrons  in  the  local  nighttime 
ionosphere.  The  sun  completes  its  journey  to  the  north  then  returns  toward  the 
magnetic  equator  [indicated  by  the  second  peak  of  the  cosine  (SM  A)  ] .  The  difference 
is  that  the  sun  has  remained  in  the  northern  magnetic  hemisphere  between  the 
peaks,  allowing  the  nighttime  ionosphere  to  become  more  ionized. 

This  difference  in  relative  sun  position  can  be  accounted  for  by  simply 
integrating  the  cosine (SMA)  over  several  past  solar  magnetic  angles,  adding  that 
value  to  the  current  SMA,  and  normalizing  to  keep  the  value  been  0  and  1 .  A  good 
fit  is  found  by  integrating  over  150  days.  The  result  of  this  integration  is  in  Figure 
38  which  shows  the  complete  modeling  of  the  circuit  probability  with  the  observed 
frequency,  in  blue,  and  the  predicted  probability  in  red. 
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Figure  38.   Probability  of  path  occurrence.   21 -day  moving  average 
circuit  probability  shown  in  blue  and  prediction  model  shown  in  red. 


Local  maxima  in  the  circuit  probability  plot  occur  when  cos(SMA)  is  zero,  i.e., 
when  the  sun  sets  on  the  magnetic  equator.  This  is  because  the  probability  of 
symmetric  anomaly  formation  is  greatest  when  solar  energy  irradiates  both  magnetic 
hemispheres  equally.  The  magnetic  hemispheres  receive  equal  solar  energy  when  the 
sun  sets  on  the  magnetic  equator. 

Local  circuit  probability  minima  occur  at  the  solstices.  As  the  setting  sun  is 
furthest  from  the  magnetic  equator  at  these  points,  asymmetric  anomaly  formation 
has  greater  probability.  Since  the  solar  position  at  the  winter  solstice  is  further  from 
the  magnetic  equator  than  at  the  summer  solstice,  the  anomaly  asymmetry  is  even 
more  pronounced,  resulting  in  a  lower  circuit  probability. 
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Absolute  symmetry  between  the  anomalies  is  seldom,  if  ever,  achieved.  It 
appears  that  there  may  be  a  symmetry  threshold,  above  which  TEP  forms  and  below 
which  it  does  not.  Symmetric  formation  is  affected  not  only  by  solar  energy,  but  also 
by  atmospheric  and  ionospheric  winds,  which  reposition  plasma.  Asymmetric 
anomaly  formation  may  account  for  non-propagation  days,  whether  the  asymmetry 
is  caused  by  solar  position  or  winds. 
D.         PATH  LOSS 

1.        Free  Space  Path  Loss 

In  line-of-sight  link  budget  calculations,  the  path  loss  is  defined  as  the  inverse 
of  the  ratio  of  received  power  to  transmitted  power,  i.e., 


Lp  = 


(Pr) 


(43) 


The  free  space  (or  spreading)  path  loss  is  a  function  of  the  transmitted  frequency  and 
the  range  from  the  transmitter  to  the  receiver, 

L  (free  space)  =  ±-±L =  -^ -,  (44) 

c2  5.7xl014 

where  f  is  in  hertz  and  R  is  in  meters. 

The  above  equation  is  in  absolute  loss,  whereas  the  path  loss  is  normally 

expressed  in  decibels  (dB).  Converting  both  sides  of  Equation  44  to  decibels, 

L^  *  32  +  20  log(0  +  20  log(R),  (45) 

where  f  is  in  MHz  and  R  is  in  kilometers. 
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2.        Received  Signal  Strength 

The  signal  power  input  to  a  receiver  is  a  function  of  the  transmitter  power,  the 
transmitter  and  receiver  antenna  gains,  and  the  path  loss.  Depending  on  the 
scenario,  there  may  or  may  not  be  control  over  power  and  gains.  The  control  over 
propagation  path  for  this  experiment  was  the  positioning  of  the  receiver  and  the 
receiving  antenna  gain. 

Because  of  the  path  variability,  the  received  TEP  signal  strength  is  not 
constant  nor  consistent.  The  total  TEP  path  losses  were  given  in  Chapter  IV  as  the 
free  space  losses,  the  coupling  losses,  and  the  reflection  losses.  The  free  space  losses 
can  be  modeled  by  assuming  a  statistical  median  distance  from  the  radio  terminals 
to  the  duct  apertures.  The  reflection  loss  is  not  known  but  presumed  to  be  very  low. 
For  modeling  purposes  it  is  assumed  to  have  a  value  of  6  dB.  As  the  radio  terminal 
positions  were  fixed  in  this  experiment,  it  was  not  possible  to  determine  coupling 
losses  resulting  from  latitudinal  and  longitudinal  terminal  variations.  However, 
using  the  known  optimum  positions,  and  the  range  of  latitudes  and  longitudes  from 
those  positions  that  allow  TEP,  coupling  factors  (inverse  of  coupling  losses)  can  be 
modeled  as  follows: 

1)  At  optimum  positions  assume  a  multiplicative  coupling  factor  of  unity 

2)  At  and  beyond  the  latitude-  and  longitude-TEP-limits  assume  a  coupling 
factor  of  zero 

3)  Assume  cosinusoidal  roll-off  at  all  points  between  1)  and  2). 
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If  it  is  assumed  that  the  distances  from  radio  terminals  to  the  duct  apertures 
have  a  median  value  of  275  km,  the  free  space  losses  can  be  computed  using  Eq.  45, 
giving  total  path  loss  as 

LT  =  LFS  +  Lc  +  LR  =  (87  +  20  log(/))  +  Lc  +  6 

(46) 
=  93+20  log(0  +  Lc  dB. 

The  free  space  path  loss  increases  with  frequency.  The  upper  limit  of  frequencies 
which  will  propagate  via  TEP  ducting  is  not  known.  The  highest  frequency  observed 
in  this  experiment,  200  MHz,  is  used  as  the  upper  limit  of  the  TEP  model.  The 
normal  minimum  HF  frequency  of  2  MHz  is  used  for  the  lower  limit. 

If  both  the  transmitter  and  receiver  are  located  at  the  optimum  magnetic 
positions  of  18°  (one  in  the  northern  hemisphere,  one  in  the  southern)  the  coupling 
losses  of  Eq.  46  would  be  zero,  but  not  at  other  latitudes.  Since  geomagnetic 
latitudes  for  TEP  are  restricted  to  14°  to  26°,  latitude  coupling  factors  are  modeled 
as  cosine  functions,  for  reasons  stated  above.  If  the  latitude  is  less  than  18°  the 
modeling  equation  for  the  coupling  factor  (CF)  is 


CFlat  -  cos  (  Wllgl 
lat  {         16 


14°  <  lat  <  18°.  (47) 


Similarly,  if  the  latitude  is  greater  than  18°,  the  equation  for  the  coupling  factor  is 

'  2n(lat-18)" 


CFla[  =  cos 


/ 


32 

The  latitude  coupling  factor  is  shown  in  Figure  39 
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18°  <  lat  <  26°.  (48) 
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Figure  39.  Latitude  coupling  factor. 
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Signal  strength  also  drops  when  the  receiver  magnetic  longitude  differs  from 
that  of  the  transmitter.  The  longitude  difference  of  this  data  set  is  6.5°  with 
substantial  reception  power,  so  a  longitudinal  limit  of  8°  is  selected.  When  signal 
strength  is  modeled  with  a  cosine  function,  to  the  limit  of  an  8°  longitudinal 
difference,  an  optimal  signal  strength  multiplicative  factor  equation  is  obtained  as 


CFlong    =   C0S 


'  2n(longV(-longRX)^ 


32 


,    0°  <  llong^-long^l  <l  8 


>RXI 


(49) 


This  factor  is  shown  in  Figure  40.  These  two  multiplicative  coupling  factors  can  be 
inverted,  expressed  in  decibels,  and  combined  as  coupling  loss  of  Equation  46  to  get 
the  total  path  loss. 
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Figure  40.  Longitudinal  coupling  factor. 


Received  signal  strength  variation  with  time  is  the  last  factor  to  be  modeled. 
The  data  shows  that  signal  strength  starts  from  zero,  builds  as  the  duct  passes  to  the 
east,  then  falls  off,  as  shown  for  five  days  in  Figure  41.  Although  not  always  true 
(e.g.,  densely  versus  sparsely  separated  ducts),  signal  strength  varies  sinusoidally 
with  time.  The  coupling  factor  equation  for  time  dependence  is  approximated  as 


CFtime 


sin 


27t(t-r0)> 

2(TC-T0)J 


(50) 


where  t  is  the  time  of  interest,  T0  is  the  circuit  onset  time,  and  Tc  is  the  circuit 
cessation  time.  Figure  42  shows  a  plot  of  the  time  coupling  factor  with  a  circuit 
onset  time  of  1830  and  a  circuit  cessation  time  of  2200.  This  multiplicative  factor 

77 


10  12 

Time  (GMT) 


Figure  41.  Signal  strength  roll-off  with  time  over  five  days.  Each  day 
is  represented  by  a  different  color. 


19.5  20  20.5  21 

Time  During  which  Circuit  is  Open 


2? 


Figure  42.   Time  coupling  factor  for  onset  at  1830  and  cessation  at 
2200  local. 
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can  likewise  be  converted  to  decibels  and  added  to  the  other  coupling  loss  factors, 
completing  the  path  loss  model  whose  modeling  equation  is 


LT  =  93  +  20  log(/)  -  10  log(CFZat  )  -  10  log(CFZar  ) 

(51) 
-  10  log(CFtap  -  10  logCCF^). 

The  prediction  of  received  signal  strength  is  found  by  subtracting  the  loss  from  the 
transmitted  power  expressed  in  decibels. 

As  validation  of  the  path  loss  model,  the  TEP  path  values  and  Eq.  51  can  be 
used  to  compute  the  expected  value  of  peak  reception  power  and  compared  with  the 
average  peak  reception  power  measurements  at  Rarotonga.  Using  the  ERP,  magnetic 
latitudes  and  longitudes,  and  transmitted  frequency  of  the  VHF  TEP  path,  the 
expected  value  of  peak  reception  power  is  found  to  be  3.51  nW.  The  peak  daily 
reception  power  values  at  Rarotonga  during  the  experiment  are  shown  in  Figure  43. 
The  mean  peak  reception  power  computed  from  the  observations  at  Rarotonga 
shown  in  this  Figure  is  3.25  nW.  Comparing  this  average  value  with  the  expected 
value  predicted  using  Eq.  51  (3.51  nW)  enhances  the  credibility  of  this  model. 

As  an  example  of  the  different  contributions  to  the  path  loss,  consider  an  onset 
time  of  1800,  cessation  time  of  2200,  transmitter  location  15 °N  magnetic,  receiver 
location  21  °S  magnetic,  and  the  difference  in  longitude  between  the  transmitter  and 
receiver  to  be  4°.  The  total  path  loss  is  computed  as 
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Figure  43.  Peak  daily  reception  power  at  Rarotonga. 
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The  result  is  graphed  in  Figure  44  for  the  duration  of  the  circuit  at  frequencies  from 
2  to  200  MHz. 
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Figure  44.  Path  loss  for  example  values. 
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VII.   THE  TRANSIONOSPHERIC  GPS  EXPERIMENT 

The  data  analyzed  in  this  chapter  are  the  GPS  signal  recordings  collected  in 
Rarotonga,  Hawaii,  and  Christmas  Island  from  March  1995  to  October  1995. 
Correlations  were  made  with  the  VHF  reception  data  collected  at  Rarotonga  over  the 
same  time  period.  The  VHF  data  are  a  continuous  monitoring  of  received  signal 
strength  from  Honolulu  and  are  considered  to  be  the  ground  truth  for  purposes  of 
this  experiment.  The  GPS  recordings  are  the  derived  total  electron  content  (TEC) 
from  all  the  visible  GPS  satellites  at  each  site,  and  for  the  positions  of  those  satellites. 

Since  the  objective  of  this  chapter  is  to  identify  and  attempt  to  predict 
scintillation  caused  by  the  equatorial  ionosphere,  the  VHF  data  must  be  used  to  help 
differentiate  between  a  scintillating  and  a  non-scintillating  GPS  signal. 
Geographically,  Christmas  Island  is  situated  virtually  beneath  the  midpoint  of  the 
propagation  path  from  Hawaii  to  the  Cook  Islands.  This  spatial  coincidence  can  be 
exploited  by  monitoring  the  GPS  signal  directly  overhead  Christmas  Island  during 
times  of  known  VHF  propagation  and  non-propagation  to  discern  the  difference 
between  the  scintillating  and  non-scintillating  GPS  signal. 
A.         TEC  OVERHEAD  CHRISTMAS  ISLAND 

If  the  Christmas  Island  data  set  is  filtered  to  exclude  all  but  GPS  signals 
arriving  from  satellites  near  directly  overhead,  these  signals  can  be  compared  with 
VHF  signals  propagating  overhead  at  the  same  time.  If  the  presence  of  ionospheric 
ducts  is  identifiable  by  using  GPS,  this  procedure  should  highlight  their  presence. 
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One  such  signal  is  that  from  GPS  satellite  2  (commonly  called  a  PRN)  on  Julian  day 
259  in  1995.  The  ground  track  of  the  signal  from  that  satellite  is  shown  in  Figure  45 


Ground  track  of  PRN  2  on  Julian  Day  259  1995 
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Figure  45.     Ground  track  showing  proximity  of  GPS  satellite  to 
Christmas  Island. 


where  the  satellite  progresses  south  to  north,  crossing  nearly  overhead  Christmas 
Island  at  2000  local  mean  time  (LMT).  With  sunset  occurring  at  approximately 
1800,  this  combination  of  satellite  position  and  time  provides  an  optimum  positional 
relationship  from  which  to  observe  the  overhead  TEC  and  the  formation  of  an 
ionospheric  duct. 

On  this  particular  day,  the  VHF  signal  from  Hawaii  to  Rarotonga  did 
propagate,  beginning  at  approximately  2000  LMT.  Recall  that  the  presence  of  the 
VHF  signal  is  indicative  of  the  existence  of  a  duct.  In  Figure  46,  the  TEC  from  PRN 
2  (shown  in  blue)  is  plotted  along  with  the  signal  strength  of  the  VHF  signal 
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TEC  of  PRN  2  and  Normalized  VHF  Signal  Strength  on  Julian  Day  259  1995 
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Figure  46.  ATEC  overhead  Christmas  Island  and  normalized  VHF 
signal  strength  at  Rarotonga.  One  TEC  unit  is  equal  to  1016 
electrons/m2. 

(normalized  to  the  maximum  TEC  value)  shown  in  red.  Notice  that  the  overhead 
crossing  time  of  2000  coincides  with  the  onset  of  the  VHF  signal  propagation,  thus 
capturing  the  TEC  response  to  the  formation  of  an  ionospheric  duct.  The  TEC  is 
seen  to  vary  between  approximately  ±0.1  TEC  units  prior  to  the  onset  of  the  VHF 
propagation.  Just  after  the  VHF  signal  begins,  the  TEC  begins  to  vary  between 
approximately  ±  1  TEC  units,  an  increase  of  1000%.  (As  noted  in  Eq.  25,  1  TEC  unit 
=  1016  electrons/m2.)  This  behavior  continues,  diminishing  as  the  VHF  signal 
strength  wanes  until  the  satellite  is  out  of  field  of  view  at  approximately  2330. 

It  is  apparent  from  the  signals  of  Figure  46  that  the  TEC  derived  from  the  GPS 
signal  is  indicative  of  the  presence  of  an  ionospheric  duct.  If  this  were  the  sole 
occurrence  of  this  type  of  TEC  response  however,  little  significance  would  be 
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attributed  to  it.  Similar  behavior  is  conspicuous  throughout  the  data  set  as  further 
evidenced  by  the  signals  shown  in  Figures  47  and  48.  Figure  47  shows  the  ground 
track  for  PRN  2  on  Julian  Day  267.  With  the  overhead  crossing  time  of  1930  the 
formation  of  another  duct  can  be  recorded  via  GPS-TEC.  This  is  illustrated  in  Figure 
48  where  the  VHF  signal  is  again  seen  to  begin  coincidentally  with  the  satellite 
passing  overhead  the  receiving  station. 

As  a  third  example,  a  day  is  chosen  on  which  there  is  VHF  propagation  but  the 
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Figure  47.  Ground  track  overhead  Christmas  Island  on  Day  267. 


satellite  overhead  position  does  not  coincide  with  the  onset  of  the  duct.  On  this  day 
the  satellite  passes  overhead  at  midnight,  2400.  Figure  49  depicts  the  onset  of  the 
VHF  signal  at  about  2230,  at  which  time  the  initiation  of  the  scintillating  TEC  signal 
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TEC  of  PRN  2  and  Normalized  VHF  Signal  Strength  on  Julian  Day  267  1995 
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Figure  48.  ATEC  and  normalized  VHF  signal  strength  on  day  267. 


TEC  of  PRN  2  and  Normalized  VHF  Signal  Strength  on  Julian  Day  201  1995 


Figure  49.  ATEC  and  normalized  VHF  signal  strength  on  day  201. 
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conforms.  This  behavior  continues  until  the  VHF  signal  ceases  at  about  0200  (shown 
as  2600  in  the  figure  for  plotting  continuity),  at  which  time  the  TEC  scintillation 
diminishes. 

For  comparison  to  the  days  on  which  there  is  no  VHF  propagation,  the  day 
after  the  previous  example  is  chosen.  On  this  day  no  VHF  signal  was  received  in 
Rarotonga  which  suggests  that  no  duct  formed  overhead  Christmas  Island.  Again 
examining  the  TEC  signal  from  PRN  2,  just  one  day  after  observing  the  response  from 
the  same  PRN  in  Figure  49,  should  yield  telling  significance  to  the  identification  of 
scintillating  TEC.  The  TEC  for  Day  202  is  shown  in  Figure  50.  The  time  of  satellite 


TEC  of  PRN  2  on  Julian  Day  202  1995 
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Figure  50.  ATEC  for  day  202.  No  VHF  signal  is  present . 

crossing  overhead  is  again  2400.   As  seen  in  the  figure,  the  amplitude  of  the  TEC 
scintillation  does  not  rise  above  the  quiescent  state  of  approximately  ±0.1  TEC  units, 
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the  apparent  noise  threshold  of  the  data.   Observation  of  the  TEC  in  its  dormant 
state,  with  no  duct  formation,  is  repeatable  throughout  the  data  set. 

The  TEC  scintillations  pictured  in  the  previous  plots  have  been  transformed 
from  slant  to  vertical  reception  angle.  The  details  of  this  transformation  will  be 
demonstrated  in  Section  E. 

It  is  apparent  from  the  foregoing  analyses  that  TEC  scintillation  is  readily 
observable  using  GPS  signals.  However,  the  methods  used  in  this  section  to  identify 
the  times  and  locations  of  scintillation  are  unacceptably  slow,  tedious,  and  labor 
intensive.  A  procedure  must  be  developed  to  more  effortlessly  identify  the  TEC 
scintillation. 
B.         FREQUENCY  DOMAIN  TEC  SCINTILLATION  ANALYSIS 

Transforming  time  series  signals  into  the  frequency  domain  for  analysis  is  a 
natural  technique  since  filters  can  be  designed  and  used  to  eliminate  or  attenuate 
unwanted  parts  of  the  signal  while  enhancing  the  desired  components.  The  initial 
step  in  analysis  in  the  frequency  domain  is  to  identify  those  frequencies  or  frequency 
bands  which  are  the  desired  aspects  of  the  signal. 

In  the  context  of  this  analysis,  it  is  desired  to  determine  the  differences 
between  the  frequency  components  of  the  quiescent  TEC  signals  and  those  of  the 
scintillated  signals.  An  obvious  way  to  compare  the  differences  is  to  transform  the 
signals  from  Figures  46,  48,  and  49,  i.e.,  the  TEC  on  those  days  during  which 
scintillation  occurs.  The  time  series  can  be  divided  into  two  parts:  the  parts  with  and 
without  scintillation.  Each  part  can  then  be  transformed  into  the  frequency  domain, 
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normalized  to  a  maximum  amplitude  of  unity,  and  compared.  It  is  anticipated  that 
the  transform  of  the  scintillating  part  will  have  frequency  domain  components  which 
are  absent  or  attenuated  in  the  non-scintillating  part  (or  vice-versa).  This  procedure 
was  employed  and  plots  of  the  results  are  shown  in  Figure  51. 

The  figure  consists  of  three  plots,  one  for  each  of  the  analysis  days,  each  of 
which  contains  the  scintillated  and  non-scintillated  spectra  for  that  day.  The 
scintillated  data  are  shown  in  red  and  the  non-scintillated  data  are  in  blue. 
Examining  each  plot,  it  is  seen  that  there  is  no  apparent  difference  between  the 
scintillated  and  the  non-scintillated  data.  This  is  true  for  any  of  the  three  days,  and 
this  result  is  found  consistently  throughout  the  data  set.  This  finding  can  be 
explained  by  considering  the  source  of  the  scintillation. 

The  TEC  scintillation  is  caused  by  the  ducts  passing  between  the  observer  and 
the  satellite.  The  ducts  themselves  consist  of  randomly  sized  tubes  imbedded  in  their 
ceilings.  Since  the  ducts  and  their  constituents  are  all  random  in  size  and  shape, 
their  effects  upon  the  signals  passing  through  them  will  be  somewhat  random.  The 
frequency  spectra  displayed  in  Figure  51  are  manifestations  of  that  randomness,  and 
are  themselves,  somewhat  random. 

This  analysis  has  shown  that  frequency  domain  techniques  for  analyzing, 
modeling,  or  identifying  the  scintillation  data  are  inadequate.  Since  the  data  are 
randomly  distributed,  the  obvious  candidates  for  further  analysis  are  statistical 
techniques. 
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Frequency  Spectra  of  Scintillated  and  non-Scintillated  TEC  on  JD  259 
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Frequency  Spectra  of  Scintillated  and  non-Scintillated  TEC  on  JD  201 
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Figure  51.  Frequency  spectra  for  three  days  containing  TEC  scintillation. 
Scintillation  spectra  shown  in  red  and  non-scintillation  spectra  shown  in  blue 
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C.         STATISTICAL  ANALYSIS  OF  TEC  SCINTILLATION 

Statistical  techniques  for  time  series  analysis  describe  how  the  data  vary  from 
different  perspectives.  Usual  descriptors  are  the  various  moments  which  can  be 
computed  from  the  data.  Moments  such  as  the  mean  and  the  variance  can  often 
differentiate  one  group  of  data  from  another.  To  help  glean  information  as  to  the 
distribution  of  the  data,  the  same  days  which  were  previously  analyzed  were  divided 
into  two  groups,  the  scintillating  and  non-scintillating.  Histograms  were  created  of 
each  to  show  their  amplitude  distributions.  These  histograms  are  shown  in  Figure 
52. 

The  top  plot  shows  the  distribution  of  the  non-scintillated  data  set  while  the 
bottom  plot  shows  that  of  the  scintillated.  Both  distributions  show  a  quasi-normal 
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Figure  52.  Histograms  of  the  scintillated  and  non-scintillated  data 
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distribution  with  a  mean  of  zero.  The  prominent  difference  between  the  two  is  the 
spread  of  the  data.  The  non-scintillated  data  are  constricted  between  approximately 
±0.25  TEC  units  while  the  scintillated  data  extend  to  approximately  ±1.75  TEC 
units.  From  these  observations  it  is  concluded  that  by  quantifying  the  spread  of  the 
data,  a  useful  scintillation  identification  tool  can  be  formulated. 

Quantifying  data  spread  is  accomplished  by  computing  the  second  and  higher- 
order  moments.  But,  since  the  data  set  is  not  stationary,  a  single  computation  of  a 
higher-order  moment  of  a  TEC  data  record,  such  as  those  shown  in  Figures  46,  48, 
and  49,  will  be  meaningless.  Instead,  it  will  be  necessary  to  compute  a  "moving 
moment",  where  a  group  of  data  is  selected,  the  moment  computed  and  stored,  then 
the  next  group  of  data  is  selected,  etc.  In  this  manner  a  new  data  record  can  be 
generated  which  represents  the  moments  of  the  data  record  with  respect  to  time. 

To  visualize  this  procedure,  the  TEC  data  from  Julian  Day  201  (Figure  49)  is 
subdivided  into  groups  of  one-half  hour  each.  The  2nd,  3rd,  and  4th  central  moments 
are  computed  for  each  half-hour  block  of  data.  The  results  are  shown  in  Figure  53 
where  the  2nd  order  moment  (variance)  is  shown  in  blue,  the  3rd  (skew)  in  green,  and 
the  4th  (kurtosis)  in  red.  The  2nd  and  4th  order  moments  appear  to  be  more 
illustrative  than  the  3rd,  and  the  2nd  appears  more  robust  than  the  4th.  Therefore,  it 
is  determined  that  the  2nd  order  central  running  moment  is  a  reliable  indicator  of  the 
presence  of  scintillation.  The  variance  has  units  of  TEC  units  squared.  The  standard 
deviation  contains  the  same  informational  value  as  the  variance,  but  its  advantage 
is  that  the  units  are  those  of  the  original  variable,  i.e.,  TEC  units.   Therefore,  the 
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Figure  53.  Variance  (blue),  Skew  (green),  and  Kurtosis  (red)  of  TEC 


scintillation  indicator  used  in  this  study  will  be  the  running  standard  deviation 
(RSD).  For  confirmation,  the  running  standard  deviation  and  the  TEC  from  which 
it  was  computed  are  plotted  together  in  Figure  54.  The  running  standard  deviation 
is  seen  to  track  the  TEC  scintillation.  The  value  of  0.25  TEC  units  for  the  running 
standard  deviation  is  subjectively  chosen  as  being  significant. 
D.         TIME-SPACE  DISTRIBUTION  OF  TEC  SCINTILLATION 

The  TEC  data  arrive  at  the  receiver  site  from  all  the  GPS  satellites  within  the 
field  of  view  of  the  receiving  antenna.  Thus,  from  one  receiving  site  at  any  given 
time,  the  TEC  scintillation  at  that  time  can  be  mapped  full-circle  around  the  site. 
This  is  achieved  by  tagging  the  received  TEC  to  the  satellite  (and  consequently  the 
azimuth  and  elevation  angle  to  it)  from  which  it  came.  Estimating  the  height  of  the 
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TEC  and  its  Running  Standard  Deviation  for  JD  201 
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Figure  54.   TEC  variation  shown  in  green  with  its  running  standard 
deviation  shown  in  blue. 

ionosphere  at  this  azimuth  and  elevation,  the  TEC  is  mapped  uniquely  to  a  three- 
dimensional  point  in  the  surrounding  sky.  Appending  the  time  stamp  to  the  three- 
dimensional  point,  a  time-space  map  of  the  TEC  can  be  constructed.  Each  point  is 
a  five-dimensional  vector:  latitude,  longitude,  altitude,  time,  and  TEC  magnitude. 
For  purposes  of  this  study  the  altitude  is  used  only  to  help  determine  the  latitude  and 
longitude.  Once  those  are  found,  the  altitude  can  be  eliminated.  But,  the  point  still 
represents  a  four-dimensional  vector. 

There  is  no  satisfactory  stationary  method  to  properly  visualize  a  collection 
of  points  in  four  dimensions.  For  example,  returning  to  Julian  Day  201,  the  RSD  can 
be  computed  for  the  two  PRNs  2  and  4.  Once  found,  the  RSD  can  be  plotted  against 
the  other  three  members  of  the  vector  one  at  a  time,  as  shown  in  Fig.  55  against  time 
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Figure  55.  Running  standard  deviation  for  PRN  2  and  PRN  4 

on  the  left  and  longitude  on  the  right  with  the  RSD  for  PRN  2  shown  in  blue  and  for 
PRN  4  in  red.  While  each  plot  displays  useful  information,  neither  shows  enough 
alone  or  together  to  gain  a  situational  awareness  of  the  environment.  The  RSD  can 
also  be  plotted  in  three  dimensions,  against  two  of  the  other  three.  The  RSD  plotted 
against  time  and  longitude  in  three  dimensions  is  shown  in  Figure  56.  Again,  PRN 
2  is  plotted  in  blue  and  PRN  4  in  red.  While  this  plot  gives  better  perspective,  the 
fourth  dimension,  latitude,  is  still  missing. 

In  order  to  properly  display  the  data  so  its  behavior  can  be  more  easily 
understood,  the  vector  dimension  must  be  reduced  from  four  to  three.  Upon 
examination  of  the  members  of  the  vector,  latitude,  longitude,  time,  and  RSD,  with 
the  additional  information  that  the  scintillation  is  a  post-sunset  event,  it  is 
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Figure  56.  Running  standard  deviation  plotted  in  three  dimensions 


ascertained  that  time  and  longitude  contain  redundant  information.  If  it  is  agreed 
that  the  parameter  that  is  actually  needed  to  describe  the  space-time  relationship  is 
time-after-sunset  (TAS),  then  time  and  longitude  can  be  combined  to  derive  that 
parameter.  This  operation  reduces  the  functional  vector  to  three  dimensions 
comprised  of  latitude,  RSD,  and  TAS. 

Before  computing  the  TAS  for  each  point  in  the  data  set,  recall  that  the 
scintillation-producing  ducts  are  aligned  along  the  magnetic  field  lines  (see  Chapter 
II).  The  magnetic  declination  angle  around  Christmas  Island  is  10.9°.  Therefore,  the 
ducts  are  not  aligned  from  true  south  to  true  north,  but  instead  are  rotated  10.9° 
clockwise.  Suppose  Christmas  Island  is  receiving  signals  from  two  PRNs  which  have 
passed  through  an  ionospheric  duct  as  shown  in  Figure  57.    There  are  two  data 
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Figure  57.  Two  GPS  signals  passing  through  different  parts  of  duct 


points  at  which  the  latitude,  longitude,  and  altitude  are  represented  by  the  impact 
points.  The  two  signals  arrive  after  being  influenced  by  the  same  scintillating  event, 
one  ionospheric  duct.  By  resolving  the  positions  of  the  scintillations,  i.e.,  latitude, 
longitude,  and  altitude,  one  finds  that  the  scintillating  cause  lies  at  similar  positive 
and  negative  latitudes,  but  two  separate  longitudes  are  returned,  211  and  214 
degrees.  While  this  is  correct  geographically,  it  does  not  reveal  that  both 
observations  are  from  the  same  event.  If  the  geomagnetic  coordinates  are  used 
instead  of  geographic,  both  signal  observations  are  properly  resolved  to  the  same 
longitude. 

With  the  data  point  positions  transformed  to  geomagnetic  coordinates,  the 
duct  shown  in  Figure  57  will  appear  as  shown  in  Figure  58.  Also  shown  in  this  figure 
is  an  arbitrarily  positioned  solar  terminator.  Computing  the  time-after-sunset,  shows 
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Figure  58.  Duct  in  geomagnetic  coordinates  showing  TAS 


the  time  from  the  upper  impact  point  to  the  solar  terminator  is  greater  than  that 
from  the  lower  impact  point.  Since  both  times  are  describing  the  same  event,  it  is 
not  desirable  to  have  two  separate  times  to  describe  the  event.  The  only  point  on  the 
solar  terminator  which  is  common  to  both  is  the  time  of  sunset  on  the  magnetic 
equator.  Accordingly,  it  is  this  time  from  which  the  TAS  is  computed,  as  shown  in 
the  figure. 

If  all  the  satellite  signals  for  an  entire  day  are  thus  transformed,  a  time-space 
map  can  be  drawn  depicting  the  TEC  scintillation  as  a  function  of  magnetic  latitude 
and  time-after-sunset  for  that  day.  This  transformation  and  mapping  was  performed 
for  four  days  using  Christmas  Island  data,  Julian  Days  187,  194,  200,  and  202.  Days 
187  and  200  experienced  VHF  propagation  while  days  194  and  202  did  not.  The 
scintillation  data  were  then  averaged  first  across  all  latitudes  and  then  across  all 
times.     These  averages  are  plotted  in  Figure  59  where  the  latitude-averaged 
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Figure  59.  Time-space  TEC  scintillation  around  Christmas  Island. 
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scintillation  is  shown  on  the  left  for  each  of  the  days  and  the  time-averaged  data  on 
the  right.  Also  on  the  left  hand  plots,  the  received  signal  strength,  in  dBm,  is  plotted 
in  green  (axis  on  right  hand  side  of  left  hand  plot).  The  VHF  signal  has  also  been 
transformed  to  time-after-sunset  on  the  magnetic  equator.  On  the  two  days  with 
VHF  propagation,  significant  scintillation  coincides  in  time.  Also,  on  the  two  days 
with  no  propagation,  scintillation  is  not  present.  On  the  right  hand  side,  the  data 
averaged  across  time  shows  the  scintillation  distribution  in  latitude.  Scintillation  is 
again  present  on  days  of  propagation  and  absent  otherwise.  The  latitudinal 
distribution  spans  from  about  -12  to  12  degrees.  It  is  confined  to  these  limits  since 
only  Christmas  Island  data  is  presented  which  limits  the  extent  of  the  field  of  view 
of  the  receiving  antenna. 

Finally,  to  gain  a  more  complete  picture  of  the  space-time  distribution  of 
scintillation,  every  GPS  signal  from  every  day  can  be  transformed  and  plotted. 
Assuming  that  every  scintillation  event  is  mirrored  in  the  opposite  hemisphere,  the 
data  from  the  southern  hemisphere  can  be  folded  on  top  of  the  northern  hemisphere 
and  vice-versa  for  visualization  enhancement.  Low-pass  filtering  in  two  dimensions 
provides  smoothing.  The  results  of  transforming,  folding,  and  filtering  the 
scintillation  data  are  shown  in  Figures  60  and  61.  The  top  plot  of  Figure  60  shows 
scintillation  as  a  function  of  time  after  sunset.  Significant  scintillation  begins  at 
approximately  1 .5  hours  after  sunset  on  the  magnetic  equator.  The  magnitude  then 
continues  to  build  until  about  5  hours  after  sunset,  at  which  time  it  diminishes  until 
about  6  hours  after  sunset  when  it  disappears.  In  the  bottom  plot,  the  scintillation 
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Figure    60.      Scintillation   standard   deviation   versus   time   and 
geomagnetic  latitude 


is  seen  to  occur  from  about  14  degrees  each  side  of  the  magnetic  equator.  Although 
there  is  still  scintillation  directly  over  the  magnetic  equator,  the  magnitude  is  not  as 
great  as  it  is  at  even  two  degrees  away.  In  the  top  plot,  surprisingly  scintillation  is 
evident  prior  to  sunset.  An  explanation  is  found  in  Figure  61,  where  the  scintillation 
before  sunset  exists  in  the  neighborhood  of  10  to  15  degrees,  the  location  of  the 
equatorial  anomalies.  This  pre-sunset  low-level  scintillation  is  caused  by  the 
anomalies  alone.  Notice  also  that  the  late-afternoon  shrinking  of  the  anomalies  is 
tracked  by  the  GPS  signals  from  -3  hours  to  -1  hour. 

On  the  other  side  of  the  solar  terminator,  the  scintillation  becomes  substantial 
about  two  hours  after  sunset.  It  spans  from  -15  to  15  degrees  with  a  maximum  at 
10  degrees.   This  continues  for  two  hours  until  four  hours  after  sunset  where  the 
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Figure  61.  Scintillation  standard  deviation  seen  from  above 

intensity  does  not  diminish  (it  increases)  but  the  latitudinal  spread  lessens  and 
withdraws  toward  the  equator.  At  five  hours  after,  the  maxima  are  found  at  about 
four  degrees.  One  final  observation  about  this  figure  is  although  meaningful 
scintillation  begins  at  about  1.5  hours  after  sunset,  there  is  no  scintillation  directly 
overhead  the  equator  for  about  another  hour. 
E.         ELEVATION  ANGLE  EFFECTS 

The  scintillation  analysis  depends  upon  the  measurements  of  TEC.  The  TEC 
is  the  integrated  sum  of  the  number  of  electrons  between  the  observer  and  the 
transmitter.  If  the  ionosphere  were  stratified,  the  fewest  number  of  electrons  on  the 
path  from  transmitter  to  receiver  would  occur  if  the  satellite  transmitter  were  directly 
overhead.  If  the  satellite  moves  from  this  position,  the  wave  travels  a  longer  path 
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through  the  ionospheric  part  of  the  journey.  A  longer  path  through  a  stratified 
ionosphere  implies  a  higher  number  of  electrons  between  the  two  sites,  a  greater 
TEC. 

This  concept  applies  to  scintillating  TEC  as  well.  The  minimum  TEC 
scintillation  will  be  measured  directly  overhead  (an  elevation  angle  of  90  degrees), 
and  will  increase  with  decreases  in  elevation  angle.  In  Figure  50,  the  TEC  of  a 
quiescent  day  was  pictured.  It  appears  to  be  uniformly  minimal  across  the  entire 
time  the  satellite  is  within  field  of  view.  However,  the  TEC  of  Figure  50  is  the 
equivalent  vertical-incidence  measurement.  In  its  unmodified  form,  the  TEC  of 
Figure  50  is  as  shown  in  Figure  62  plotted  against  the  elevation  angle  at  which  it  was 
collected.  Here  it  is  seen  that  although  the  TEC  is  from  an  undisturbed  ionosphere, 
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Figure  62.  Variation  of  TEC  as  function  of  elevation  angle. 
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the  TEC  increases  as  the  elevation  angle  decreases  from  its  most-near  overhead 
position,  with  radical  changes  as  the  elevation  angle  approaches  twenty  degrees. 
This  same  elevation  angle  dependence  is  ubiquitous  throughout  the  data  set. 

The  scintillation  that  has  been  analyzed  and  discussed  is  relevant  only  to  the 
overhead  (elevation  angle  90  degrees)  position.  At  other  elevation  angles  the  TEC 
and  its  scintillation  values  will  be  increased  based  upon  elevation  angle.  To  allow 
for  that  effect,  a  compensation  factor  must  be  included  in  any  scintillation 
computation  algorithm.  A  reasonable  compensation  factor  is 


TECvertical    "   ^Cslant 


1  -  e 


A  -  18.5 
10 


(53) 


(where  A  is  the  elevation  angle  in  degrees),  and  is  plotted  in  Figure  63.  This  factor 
was  used  to  transform  the  TEC  data  which  have  been  used  in  the  foregoing  analyses. 
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Figure  63.  Slant  TEC  to  Vertical  multiplicative  compensation  factor. 
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When  the  receiver  is  aboard  a  space-borne  vehicle,  the  elevation  angle  can  be 
estimated  using  a  topocentric  system  at  the  estimated  location  of  the  transmitter. 
This  estimated  elevation  angle  can  then  be  used  in  Eq.  53. 
F.         TEC  SCINTILLATION  EFFECTS  UPON  GEOLOCATION 

In  Section  D  it  was  determined  where  and  when  scintillation  will  occur. 
However,  its  effects  at  those  locations  were  not  addressed.  The  scintillation  can 
cause  deviations  in  time,  frequency,  and  angle  measurements.  Using  the  derivations 
of  Chapter  IV  the  errors  caused  by  each  of  those  factors  will  be  discussed. 

1.        Time  Difference  of  Arrival  Errors 

In  Equation  30,  it  was  determined  that  the  time  of  arrival  of  the  received 
signal  is  altered  by  the  TEC  using  the  relationship 

Af  =   1.34*10"  .  rec  m  1.34*10>  .  rac  UnUs 
f2  f2 

Time  difference  of  arrival  geolocation  systems  resolve  distance  from  each  receiver  to 
the  transmitter  by  multiplying  time  by  the  speed  of  light.  Multiplying  Eq.  54  by  c 
results  in  the  distance  error  of  cAt.  A  percentage  error  in  the  distance  measurement 
can  be  computed  based  upon  the  height  of  the  space  platform.  For  two  platforms, 
one  at  1000  km  height  and  the  other  at  20000  km  height,  the  percentage  error  at  one 
TEC  unit  is  computed  and  shown  in  Figure  64.  The  percentage  error  is  greater  for 
the  space  platform  closer  to  earth.  For  example,  assume  the  frequency  of  the 
transmitted  frequency  is  10  MHz,  the  receiving  platform  is  directly  overhead  at  1000 
km,  with  one  TEC  unit  of  scintillation.  From  the  plot,  the  percent  error  in  distance 
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Figure  64.  Distance  measurement  errors  from  two  satellites 


measurement  will  be  approximately  one  percent.  In  Figure  60,  the  TEC  scintillation 
is  seen  to  be  approximately  three  TEC  units  for  most  of  the  post-sunset,  pre-midnight 
sector.  In  this  case,  the  example  cited  above  would  increase  to  a  three  percent  error 
in  distance  measurement.  If  the  elevation  angle  from  the  transmitter  were  20 
degrees,  Eq.  53  can  be  used  to  find  the  scintillation  at  this  elevation  angle,  i.e., 


TEC 


TEC 


vertical 


slant 


1 


A        18.5 
10 


(55) 


With  TEC^rtfcd  at  three  TEC  units  and  A  at  twenty  degrees,  the  TECslant  is  found  to  be 
approximately  18  TEC  units,  with  a  resulting  percentage  error  of  approximately  50 
percent. 
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By  this  analysis,  the  optimum  observation  position  is  directly  overhead.  As 
the  position  of  the  observer  moves,  the  errors  in  the  resolution  progressively  worsen. 
The  amount  of  degradation  is  algorithm  specific.  Time  and  distance  errors  predicted 
here  must  be  applied  to  the  user's  algorithm  to  determine  individual  effects. 

2.         Frequency  Difference  of  Arrival  Systems 

From  Equation  32,  the  change  in  frequency  resulting  from  the  change  in  TEC 
is  given  by 


1  34rl09    d 
A/  =    1"ytyiu    ±  (JEC  Units)  Hz. 

f  dt 


(56) 


As  an  example  of  change  in  TEC  with  respect  to  time,  the  TEC  of  PRN2  on  Julian  Day 
259  shown  in  Figure  46  is  expanded  around  time  2100  as  shown  in  Figure  65.  The 
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Figure  65.  Change  in  TEC  with  change  in  time 
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plot  shows  a  change  in  TEC  of  about  3  TEC  units  during  a  time  change  of 
approximately  0.14  hours,  or  0.006  TEC  units  per  second.  Substituting  this  value  for 
d/dt(TEC  Units)  in  Eq.  56,  the  change  in  frequency  for  this  case  can  be  computed. 
Figure  66  shows  this  frequency  deviation  as  a  function  of  transmitted  frequency.  For 
a  transmitted  frequency  of  10  MHz,  the  change  in  frequency  at  the  receiver  is  seen 


Frequency  (hertz) 

Figure  66.  Frequency  deviation  caused  by  scintillation 

to  be  about  one  hertz.  Based  on  the  previous  derivation  based  upon  three  TEC  units, 
at  20  degrees  elevation,  the  frequency  change  is  about  18  hertz. 

These  values  must  be  tested  in  the  user's  algorithm  to  determine  the  extent 
of  the  degradation  to  geolocation. 
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3.        Angle  of  Arrival  Systems 

To  use  wedge  refraction  theory,  the  change  in  TEC  with  respect  to  horizontal 
distance  is  required.  The  time  derivative  of  Figure  65  can  be  used  to  estimate  this 
spatial  derivative.  The  TEC  time  derivative  was  found  to  be  approximately  3  TEC 
units  in  0.14  hours.  On  the  equator,  one  hour  equates  to  15  degrees  of  longitude 
which  is  equivalent  to  1665  kilometers.  The  spatial  derivative  is  estimated  as  3  TEC 
units  per  233  kilometers  or  13  /zTEC  units  per  meter.  From  Eq.  33 


40.3xl016    d    ,^„  TT  .  x 
a  =  —  (TEC  Units), 

f2  dx 


(57) 


where  a  is  in  radians.  Substituting  for  the  spatial  derivative,  the  angle  deviations  can 
be  computed  as  a  function  of  frequency  as  shown  in  Figure  67.   At  10  MHz,  the 
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Figure  67.  Angular  deviation  due  to  spatial  TEC  variation 
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angular  change  is  0.1  radian,  or  5.7  degrees.  At  1000  km  altitude,  this  equates  to 
about  25  kilometers  on  the  ground.  As  before,  this  plot  is  for  an  overhead  receiver. 
Receivers  at  positions  other  than  the  optimum  suffer  even  more  severe  errors. 
G.         SCINTILLATION  PREDICTION 

The  prediction  of  the  ionospheric  TEC  scintillation  can  be  divided  into  two 
different  areas.  First,  the  prediction  of  the  probability  that  scintillation  will  exist,  and 
second,  when  scintillation  is  present,  where  it  will  occur  and  with  what  magnitude. 

1.  Probability  of  Scintillation 

Analysis  of  the  data  in  Figures  46,  48,  49,  and  50  shows  that  the  presence  of 
the  VHF  signal  implies  scintillation  and  vice-versa.  This  is  shown  for  a  few  examples. 
Thus,  a  prediction  of  one  is  a  prediction  of  the  other.  Prediction  of  VHF  signal 
propagation  was  developed  in  Chapter  VI.  The  MATLAB  program  used  to  compute 
the  VHF  probability  has  been  modified  slightly  to  predict  scintillation  and  is  included 
as  Appendix  E. 

2.  Space-Time  Distribution 

The  distribution  prediction  modeling  is  based  upon  the  results  shown  in 
Figures  60  and  61.  The  top  plot  of  Figure  60  shows  the  magnitude  of  the  scintillation 
based  upon  time  after  sunset  can  be  modeled  with  two  exponentials,  one  to  model 
the  rise  and  peak  and  the  other  to  model  the  fall.  The  equations  of  these  two 
exponentials  are 

/  TAS  -  \\ 

timemagl  =  3.5  •  { 1  -  e      09    j ,         1  <  TAS  $  5.5, 
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and  zero  otherwise  for  the  first,  and 


timemag2  =  3.5  •  \e 


TAS  -  5.5  ^ 
0.3 


5.5  <  TAS  <  7, 


(59) 


and  zero  otherwise  for  the  second  where  TAS  is  time  after  sunset.  Scintillation  with 
respect  to  time  after  sunset  is  the  sum  of  timemagl  and  timemag2.  A  plot  showing 
both  the  top  plot  of  Fig.  60  and  this  sum  is  shown  in  Figure  68.  Note  that  since  this 
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Figure  68.  Time  after  sunset  observation  and  model. 


is  a  nighttime  effects  model,  all  effects  prior  to  one  hour  after  sunset  have  been 
eliminated.  The  purpose  of  this  mathematical  model  is  to  predict  the  scintillation 
versus  time.  Scintillation  is  also  a  function  of  latitude,  as  shown  in  Figures  60  and 
61. 
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To  derive  mathematical  descriptions  of  the  scintillation  variation  with  respect 
to  latitude,  the  bottom  plot  of  Figure  60  and  Figure  61  will  be  used.  First,  the  peak 
latitude  as  a  function  of  time  will  be  patterned.  The  greatest  latitudinal  extent  is 
about  15  degrees  occurring  approximately  two  hours  after  sunset.  The  latitude  peak 
for  this  time  is  12  degrees.  The  peak  latitude  decreases  from  there  until  six  hours 
after  sunset.  A  cosine  fit  to  this  peak  and  fall-off  is 


peaklat  =  12  cos 


'  2n(TAS  -  2)) 


17.5 


1  <;  TAS  <  7. 


(60) 


From  the  latitude  of  peak  scintillation  at  any  given  time,  the  magnitude  falls 
off  steeply  at  latitudes  greater  than  peak,  and  gradually  at  latitudes  less.  For  the 
greater  latitudes,  the  magnitude  falls  as  a  cosine  with  zero  magnitude  attained  at  six 
degrees  above  peak  latitude,  i.e., 


abovemag  =  0.5 


/ 

/ 

1 

+ 

cos 

V 

I 

2n(peaklat  -  maglat)^ 
12 


(61) 


peaklat  ^  maglat  £  peaklat  +  6,  where  maglat  represents  geomagnetic  latitudes,  in 
this  case  greater  than  the  peak  latitude.  The  magnitude  at  latitudes  less  than  peak 


is 


belowmag  =  0.5 


/ 

/' 

1 

+ 

COS 

V 

\ 

2iz(peaklat  -  maglat) 
32 


\  \ 


// 


(62) 


0  <>  maglat  <  peaklat. 

The  predicted  scintillation  magnitude  at  a  given  TAS  and  geomagnetic  latitude 
can  be  determined  by  the  following  steps: 
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1)  Determine  the  magnitude  with  respect  to  time  using  either  Eq.  58  or  59. 
Call  this  value  timemag. 

2)  Use  Eq.  60  to  determine  the  location  of  peak  latitude,  called  peaklat. 

3)  Either  Eq.  61  or  62  is  used  to  determine  the  modifying  factor  based  upon 
latitude.  Call  this  value  latfactor. 

4)  The  predicted  value  of  scintillation  is  then  the  product  of  timemag  and 
latfactor. 

A  MATLAB  program  to  compute  these  values  in  contained  in  Appendix  F.  Figure  69 

shows  a  plot  with  all  relevant  times  and  latitudes  and  Figure  70  shows  an  overhead 


view. 


Geomagnetic  Latitude  (degrees) 


■20        1 
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Figure  69.  Plot  of  predicted  TEC  scintillation 

The  equations  and  plots  above  are  for  overhead  observer  only.  If  the  elevation 
angle  from  earth  to  satellite  is  not  90  degrees,  then  the  results  obtained  above  must 
be  modified  from  vertical  to  slant  TEC  using  Equation  55. 
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Figure  70.  Scintillation  prediction  from  above. 

The  magnitudes  of  the  scintillation  predictions  are  based  upon  a  data  set 
collected  at  sunspot  minimum.  The  scintillation  dependence  upon  sunspot  number 
cannot  be  determined  from  this  data  set,  but  it  is  anticipated  that  the  observations 
here  are  the  "best  case"  conditions. 
H.         NON-TEP  SCINTILLATION 

It  has  been  shown  that  TEP  and  low-latitude  scintillation  are  separate, 
independent,  but  correlated  manifestations  arising  from  the  Rayleigh-Taylor 
instability  of  the  anomalies.  When  TEP  occurs,  an  ionospheric  duct  is  the  conduit 
which  permits  the  propagation.  Due  to  the  irregularities  in  the  duct  ceiling, 
scintillation  of  a  transionospheric  signal  coexists  with  the  TEP. 
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There  exists  the  possibility  of  the  Rayleigh-Taylor  instability  developing  within 
each  anomaly  without  formation  of  ionospheric  ducts.  Bubbles  could  rise  within  the 
anomalies  but  never  elongate  sufficiently  to  adjoin.  Alternatively,  a  bubble  could 
form  in  one  anomaly  but  not  the  other.  In  these  circumstances  no  TEP  would 
transpire,  but  scintillation  would. 

The  measure  of  scintillation  magnitude  adopted  here  has  been  the  running 
standard  deviation  of  TEC.  The  peak  RSD  for  each  day  of  the  experiment  is  shown 
in  Figure  71.  The  average  peak  RSD  for  the  entire  record  is  0.99  TEC  units. 
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Figure  71.   Peak  RSD  for  each  day  of  the  GPS  experiment.   The  line 
between  days  218  and  230  indicates  an  absence  of  data  on  those  days. 
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The  peak  RSD  for  days  with  and  without  TEP  can  be  plotted  separately  to 
determine  if  scintillation  is  significant  on  days  without  TEP.  Figure  72  shows  the 
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Figure  72.  Peak  RSD.  The  top  plot  shows  the  RSD  for  days  when  TEP 
occurred,  while  the  bottom  depicts  the  values  when  TEP  did  not  occur. 


peak  RSD  over  the  record  divided  into  days  when  TEP  occurred,  shown  in  the  top 
plot,  and  days  when  it  did  not,  on  the  bottom.  The  RSD  on  TEP  days  shows  a 
maximum  of  about  seven  TEC  units  while  the  non-TEP  plot  has  a  maximum  of  about 
0.7.  The  average  peak  RSD  for  TEP  days  is  1.27,  and  for  non-TEP  0.29  TEC  units. 
While  it  cannot  be  stated  that  scintillation  does  not  exist  on  days  when  TEP 
is  not  present,  it  is  apparent  that  its  severity  is  greatly  reduced.  The  average  value 
of  0.29  TEC  units  on  non-TEP  days  barely  exceeds  the  significance  threshold 


117 


established  earlier  of  0.25  TEC  units.  It  is  therefore  concluded  that  scintillation  does 
not  occur,  in  a  practical  sense,  without  the  formation  of  a  TEP  conducting 
ionospheric  duct. 
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VIII.  THE  ANCILLARY  EXPERIMENTS 

Two  ancillary  investigations  were  performed  which  supplement  the  VHF 
transequatorial  and  GPS  transionospheric  experiments.  The  first  supplemental 
experiment  was  the  measurement  of  transequatorial  propagation  in  the  opposite 
direction  to  the  VHF  propagation.  The  other  was  the  evaluation  of  airglow  overhead 
Christmas  Island  which  was  taken  to  augment  the  GPS  scintillation  data  there. 
A.         THE  IONOSONDE  EXPERIMENT 

An  oblique  incidence  digital  ionosonde  was  activated  from  Rarotonga  to  Oahu 
over  the  time  period  26  August  to  25  September  1995.  The  purpose  of  the 
experiment  was  two-fold:  (1)  determine  the  bi-directionality  of  the  VHF  TEP,  i.e., 
determine  if  propagation  in  one  direction  implies  propagation  in  the  other  and  (2) 
allow  observation  of  afternoon-TEP. 

The  information  supplied  by  an  ionogram  consists  of  a  plot  of  the  time-of- 
flight  for  the  frequencies  which  propagate.  In  order  to  make  some  sense  of  the  time- 
of-flight  numbers,  an  analysis  of  the  geometry  of  the  path  must  be  undertaken. 
Using  the  geometrical  distances  of  different  flight  profiles  and  a  propagation  velocity 
of  c,  estimated  times-of-flight  can  be  obtained. 

The  two  forms  of  TEP,  afternoon  and  evening,  arrive  at  the  receiver  from 
different  paths.  The  afternoon-TEP  path  allows  reflection  from  one  anomaly  to  the 
other,  then  to  the  receiver.  In  order  to  compute  the  path  length  and  propagation 
time,  both  the  heights  and  positions  of  the  anomalies  must  be  estimated.    The 
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geometry  of  the  afternoon-TEP  path  is  pictured  in  Figure  73.  The  path  is  divided  into 
three  segments:  Px  from  transmitter  to  first  reflection  point,  P2  from  the  first  to  the 


Figure  73.  Afternoon-TEP  path  geometry, 
second  reflection  point,  and  P3  from  second  reflection  point  to  receiver.  Simple 
trigonometry  using  the  radius  of  the  earth  (Re),  hlf  h2,  and  the  angles  Qlt  02,  and  63 
can  be  used  to  solve  for  the  three  path  lengths.  Since  the  anomaly  positions  are  not 
known  a  priori,  the  values  of  6X,  02,  03,  hlt  and  h2  must  be  estimated.  For  typical 
values  of  anomalies  positioned  at  ±  12°  geomagnetic  latitude,  62  is  24°.  The  angles 
0X  and  03  can  be  found  using  the  transmitter  and  receiver  latitudes.  Reflection 
heights  can  be  estimated  as  h:  =  h2  =  250  km.  Using  these  values,  a  path  length  of 
approximately  4900  km  is  found  with  a  time-of-flight  of  16.3  milliseconds.  This 
represents  the  flight  time  for  the  low  ray  reflected  from  the  anomalies.  There  exists 
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the  possibility  of  a  high  ray,  as  well.  The  time-of-flight  of  the  high  ray  will  be  greater 
than  that  of  the  low  ray,  except  at  the  maximum  propagation  frequency,  where  they 
merge. 

For  evening  TEP,  if  the  ray  path  follows  a  magnetic  field  line  from  transmitter 
to  receiver,  an  estimate  of  the  path  length  can  be  attained.  The  length  of  the 
magnetic  field  line,  from  Equation  1,  produces  a  path  length  of  approximately  5025 
km  with  a  flight  time  of  16.7  ms. 

Using  these  estimates  of  time-of-flight  over  paths  of  the  two  types  of  TEP, 
selected  ionograms  recorded  during  the  experiment  were  examined.  The  following 
ionograms  were  selected  from  the  recordings  made  on  September  25.  The  first  one, 
taken  at  local  time  1527,  is  a  normal,  HF-limited  ionogram,  shown  in  Figure  74.  The 
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Figure  74.  HF  ionogram  at  1527. 
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maximum  observed  frequency  (MOF)  is  about  30  MHz.  The  time-of-flight  is 
consistent  with  that  of  high  altitude  propagation  with  possibilities  of  multiple  hops 
evident.  This  is  not  TEP,  but  indicates  a  high  density  ionosphere,  reflecting  HF. 

In  the  second  ionogram,  Figure  75,  taken  10  minutes  later,  the  same  HF 
reflected  traces  are  still  visible,  but  a  new  trace  reflecting  much  higher  frequencies 
appears.    This  new  trace,  arriving  at  approximately  16.3  ms,  is  the  emerging 
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afternoon-TEP  with  a  MOF  of  about  40  MHz.  The  time-of-flight  is  consistent  with  the 
theoretical  time-of-flight  for  afternoon-TEP  derived  above.  Both  the  low  and  high 
rays  are  seen  on  the  trace.  The  TEP  trace  continues  to  strengthen  and  allows  higher 
propagation  frequencies,  shown  in  Figure  76.  The  MOF  is  about  46  MHz  at  1632 
local  time. 
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Figure  76.  Afternoon-TEP  extending  to  46  MHz. 
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The  afternoon-TEP  trace  continued  until  after  sunset.  At  1852  local  time,  the 
recorded  ionograms  begin  a  radical  alteration  as  shown  in  Figure  77.  The  afternoon- 
TEP  trace  has  all  but  disappeared  but  the  trace  at  approximately  16.7  ms  is  extending 
into  the  VHF  region.  Twenty  minutes  later,  the  continued  extension  of  this  trace 
approaches  50  MHz,  as  shown  in  Figure  78.  The  time-of-flight  of  the  frequencies 
observed  on  this  trace  is  also  consistent  with  that  derived  for  evening-TEP  and  shows 
the  emergence  of  evening-TEP.  This  behavior  continues  until  the  layer  reflections 
have  all  but  disappeared  and  the  MOF  continues  to  increase,  as  shown  in  Figure  79. 
The  reception  of  evening-TEP  continues  throughout  the  evening,  but  in  Figure  80, 
near  midnight,  signals  begin  to  weaken,  as  the  propagation  mechanism  is  dissipating. 
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Figure  77.  Afternoon-TEP  trace  disappears  at  1852. 
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Local  time  1922 
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Figure  79.  Evening-TEP. 
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Figure  80.  Dissipating  evening-TEP. 
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An  appreciation  for  die  frequencies  the  two  TEP  types  will  support  is  obtained 
from  the  maximum  observed  frequency  of  each  ionogram,  shown  in  Figure  81.  The 
maximum  observed  frequency  of  evening-TEP  was  65  MHz  while  the  afternoon-TEP 
maximum  was  about  45  MHz. 
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Figure  81.  Maximum  observed  frequencies  over  24  hours  at  Oahu. 


The  bi-directionality  of  the  path  is  established  if  for  each  night  of  VHF 
evening-TEP  there  is  also  a  TEP  on  the  sounder  circuit.  The  onset  and  cessation 
times  for  both  paths  for  the  days  when  the  sounder  was  operational  are  plotted  in 
Figures  82  and  83.  The  plots  show  that  both  the  start  and  stop  times  are 
comparable.  More  important,  the  paths  are  bi-directional  with  the  exceptions  of 
Julian  days  239  (south-to-north  only)  and  260  (north-to-south  only). 
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Figure  82.  VHF  and  sounder  onset  times. 
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Figure  83.  VHF  and  sounder  cessation  times. 
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B.         THE  AIRGLOW  EXPERIMENT 

An  airglow  monitor  was  placed  at  Christmas  Island  from  14  September 
through  2  October  1995.  The  purpose  of  observing  the  airglow  variations  overhead 
Christmas  Island  was  two-fold.  The  first  and  primary  objective  was  to  obtain  an 
independent  verifier  of  the  GPS-TEC  scintillation  data.  The  second  objective  was  to 
appraise  airglow  monitoring  as  a  real-time  evaluator  of  the  presence  of  ionospheric 
ducts. 

Ionospheric  scintillation  was  shown  in  Chapter  VII  to  be  observable  using  GPS. 
Simultaneous  observations  of  the  ionosphere  directly  overhead  Christmas  Island  via 
GPS  signal  scintillation  and  airglow  variation  provides  independent  and  useful 
verification  of  the  findings  of  Chapter  VII. 

Comparison  of  the  two  data  sets  is  not  trivial.  Since  the  airglow  data  consists 
of  zenith  intensity  measurements,  the  GPS  data  to  correlate  it  with  must  come  from 
directly  overhead,  also.  This  limits  the  amount  of  data  available  for  comparison  as 
many  GPS  satellites  do  not  cross  overhead  on  a  given  day.  Furthermore,  satellites 
which  do  cross  overhead  pass  quickly,  making  comparison  even  more  difficult. 

GPS  signals  observed  solely  at  zenith  do  not  provide  enough  information  for 
comparison.  The  satellite  signal  must  be  observed  for  some  period  of  time  before 
and  after  the  time  of  overhead  crossing.  The  filtering  criterion  used  in  this  analysis 
for  GPS  signal  comparison  was  ±4°  from  zenith  in  both  latitude  and  longitude. 
While  this  window  is  narrow  enough  to  ensure  spatial  correlation,  observations  off 
zenith  must  be  translated  in  time  to  transform  them  to  the  overhead  position.  The 
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airglow  and  GPS  signals  are  not  calibrated  to  the  same  units  so  absolute  magnitude 
comparisons  are  not  useful.  Instead,  both  signals  are  normalized  so  that  relative 
magnitude  changes  can  be  compared.  After  transformation  and  normalization,  the 
airglow  and  GPS  signals  can  be  plotted  together  for  correlation. 

The  plot  of  airglow  and  GPS  data  from  Julian  Day  271  in  Figure  84  show 
airglow  intensity  variations  in  blue  while  the  signals  from  several  GPS  satellites  are 
in  green.  On  this  night,  TEP  did  occur  and  scintillation  is  evident  on  the  GPS  signals. 
Note  the  significant  variability  in  the  airglow  intensity  signal.  Although  the  airglow 
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Figure  84.  Normalized  Airglow  and  TEC  during  TEP. 

and  GPS  signals  do  not  closely  track  each  other,  it  is  evident  that  similar  processes 
created  both  and  that  they  are  correlated. 
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Figure  85  shows  a  similar  plot  for  Julian  Day  257  when  TEP  did  not  occur. 
Only  one  satellite  passed  overhead  during  the  time  of  airglow  data  recording  but  its 
waveform  in  green  follows  the  airglow  trace  in  blue.  Equally  as  important  is  the  fact 
that  the  structure  formed  during  TEP  (shown  in  the  measurements  of  Figure  84),  is 
absent  in  this  figure,  when  there  is  no  TEP. 
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Figure  85.  Normalized  Airglow  and  TEC  with  no  TEP. 


It  is  evident  from  these  figures  that  airglow  intensity  is  correlated  with  GPS 
TEC  measurements.  It  is  also  seen  that  airglow  observations  are  a  reliable  indicator 
of  the  presence  of  ionospheric  ducts.  However,  complete  correlation  of  the  two 
signals  was  impossible  because  only  the  zenith  airglow  data  were  analyzed.  This 
limitation  eliminated  most  of  the  GPS  measurements  from  the  correlation  data  pool, 
permitting  only  a  cursory  comparison. 
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IX.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  two  major  objectives  of  this  research  were  to  attain  an  ability  to  predict 
(1)  the  when  and  where  of  transequatorial  propagation  and  (2)  the  temporal  and 
spatial  distribution  of  low-latitude  transionospheric  signal  scintillation.  These  goals 
were  met  using  radio  circuits,  GPS  receivers,  airglow  monitoring,  and  an  oblique 
ionosonde.  The  primary  contributing  factor  to  the  predictive  algorithms  was  found 
to  be  the  relative  solar  intensity. 
A.         TRANSEQUATORIAL  PROPAGATION  PREDICTION 

From  observations  nd  analysis  of  transequatorial  VHF  propagation,  a  list  of 
the  physical  parameters  which  exercise  control  over  propagation  variables  has  been 
identified.  Propagation  variables  which  were  identified,  and  for  which  predictions 
were  produced,  are  the  circuit  onset  time,  the  circuit  cessation  time,  the  probability 
of  the  existence  of  the  circuit,  and  the  path  loss.  The  physical  parameters,  used  as 
predictors,  controlling  these  variables  are  the  10.7  cm  solar  flux,  the  magnetic 
declination  angle,  and  the  solar  magnetic  angle,  all  related  to  the  intensity  and 
relative  position  of  the  sun. 

Using  algorithms  derived  in  Chapter  VI,  a  user  must  simply  supply  the  latitude 
and  longitude  of  a  transmitter,  the  frequency,  the  solar  flux  value,  and  the  date.  For 
a  given  receiver  location,  the  outputs  are: 

1)  the  probability  of  a  TEP  circuit  developing, 

2)  the  predicted  starting  and  stopping  times  given  that  the  circuit  exists,  and 
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3)  the  path  loss. 

The  propagation  path  was  shown  to  be  bi-directional  over  28  of  30  days,  93% 
of  the  time.  On  the  two  days  when  bi-directionality  did  not  exist,  one  day  allowed 
propagation  north-to-south  and  the  other  the  south-to-north,  eliminating  different 
circuit  gain  margins  as  a  cause.  While  greater  than  90%  is  within  acceptable  limits 
of  ionospheric  predictability,  it  is  not  known  why  signals  will  propagate  in  one 
direction  and  not  in  the  other.  Possibly  the  magnetic  field  renders  the  ducts  non- 
reciprocal  so  that  bi-directionality  was  present  on  all  the  days  but  signal  strength  in 
one  direction  was  severely  attenuated.  The  science  of  low-latitude  ionospheric  ducts 
is  still  maturing. 
B.         SCINTILLATION  AND  GEOLOCATION 

Global  positioning  receivers  have  proven  to  be  an  invaluable  asset  in 
monitoring  both  the  quiescent  and  disturbed  ionosphere.  By  capitalizing  on  the 
coherent  frequencies  of  the  GPS  signal,  the  traveling  disturbances  in  the  equatorial 
region  can  be  observed,  quantified,  and  predicted. 

The  occurrence  of  scintillation  was  found  to  be  very  highly  correlated  with 
TEP.  No  cases  of  severe  running  standard  deviation  scintillation  were  observed  in 
the  absence  of  TEP. 

The  predictability  of  low-latitude  scintillation,  both  in  time  and  in  space,  has 
been  demonstrated.  Mathematical  relationships  have  been  developed  to  help  predict 
the  scintillation  value  at  an  arbitrary  point  in  time  and  location  in  the  low  latitude 
ionosphere. 
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The  impact  of  low-latitude  ionospheric  scintillation  upon  the  geolocation 
solution  was  found  to  vary  from  catastrophic  to  minimal,  depending  on  time, 
frequency,  and  radio  signal  penetration  angle  through  the  ionosphere.  Highlights  of 
the  findings  are  summarized  below. 

Geolocation  at  low  latitudes  can  be  severely  affected  by  the  nighttime 
ionosphere.  The  low  latitude  scintillation  extends  from  about  -20  to  20  degrees 
geomagnetic  latitude.  The  effect  on  the  geolocation  solution  is  three-fold:  temporal, 
frequency,  and  elevation  angle. 

1.  Temporal  Impact 

In  this  report  only  nighttime  ionospheric  effects  were  examined.  Nighttime 
effects  are  the  most  severe  from  one  to  seven  hours  after  sunset  at  the  magnetic 
equator.  The  recurrence  rate  of  scintillation  is  a  function  of  sun  position. 

2.  Frequency  Dependence 

The  effect  of  frequency  on  a  geolocation  solution  depends  on  which  system 
(TDOA,  AOA,  or  FDOA)  is  employed.  In  general,  when  scintillation  is  present, 
signals  below  one  megahertz  are  so  severely  disturbed  that  geolocation  using  any  of 
the  methods  is  impossible,  where  the  signal  path  is  through  the  low-latitude 
ionosphere.  For  microwave  frequencies  (above  10  gigahertz),  minimal  disturbance 
occurs  for  any  of  the  geolocation  methods.  Frequencies  between  these  two  limits  are 
affected  by  varying  levels  of  severity. 

The  preferred  position  for  the  geolocation  receiver  is  directly  above  the 
emitter,  at  a  90  degree  elevation  angle  from  the  terrestrial-based  transmitter.  For 
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comparison,  an  elevation  angle  of  90  degrees  is  assumed,  with  the  receiver  at  1000 
km  altitude.  At  three  TEC  units  of  scintillation  standard  deviation,  errors  generated 
from  the  scintillated  portion  of  the  geolocation  solution  are: 

Overhead  Geolocation  Errors 


Frequencv  rMHz^ 

TDOA  rkm) 

FDOA  THz^) 

AOA  fknri 

1 

1200 

24 

7860 

10 

12 

2.4 

79 

100 

0.12 

0.24 

0.8 

1000  0.0012  0.024  0.008 

3.        Elevation  Angle  Dependence 

The  elevation  angle  is  defined  from  the  point  of  view  of  the  terrestrial 
transmitter,  with  the  vertical  transmission  direction  at  90  degrees  and  a  horizontal 
launch  at  zero-degree  elevation  angle.  The  degree  of  degradation  of  the  arriving 
signal  is  highly  dependent  upon  the  elevation  angle.  A  receiver  positioned  directly 
above  the  transmitter  has  the  maximum  resolution,  with  lower  elevation  angles 
producing  worse  solutions.  For  example,  at  40  degrees  elevation,  the  scintillation 
magnitude  increases  by  13%,  at  30  degrees  by  50%,  and  at  20  degrees  by  700%. 
This  elevation  angle  effect  is  illustrated  in  the  computations  of  Paragraph  2,  repeated 
below  for  20  degrees  elevation. 
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Geolocation  Errors  at  Twenty  Degree  Elevation  Angle 
Frequency  fMHz)  TDOA  fkm^)  FDQA  (Hz)  AOA  (km) 

1  8680  174  56000 

10  87  17.4  560 

100  0.87  1.74  6 

1000  0.009  0.174  0.06 

C.         EXTENSIONS  OF  CURRENT  STUDY 

For  extending  the  knowledge  and  predictive  ability  gained  from  this  study, 
several  additional  investigations  are  needed.  Seven  of  these  follow. 

1.  Longitudinal  Extension 

A  useful  supplement  to  the  data  analysis  of  this  study  would  be  to  analyze 
existing  GPS  data  at  other  longitudes  in  the  Pacific  during  the  same  time  period  as 
that  of  this  data  set.  Data  sets  that  could  be  helpful  are  those  from  Tahiti, 
Central/South  America,  and  Kwajalein  Island.  Analysis  could  assist  in  determining 
longitudinal  correlation  of  nighttime  equatorial  scintillation.  If  fruitful,  an 
observation  of  scintillation  in  Central  America  could  provide  information  on  the 
probability  of  scintillation  near  Hawaii  a  few  hours  later.  Similarly,  scintillation 
observation  in  Christmas  Island  might  help  predict  scintillations  over  the  Phillipines 
on  the  same  evening. 

2.  Scintillation  Monitoring  at  Higher  Sampling  Rate 

No  frequency  distinctions  were  found  in  the  scintillation  data.  One  limitation 
with  the  data  set  is  the  sampling  period  at  Christmas  Island  and  Rarotonga  was  six 
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seconds.  The  JPL  data  is  even  more  limited,  with  sampling  at  30  seconds.  A  six- 
second  sampling  period  restricts  the  highest  observable  frequency  to  80  millihertz 
which  may  be  too  low.  One  solution  is  to  collect  additional  data  at  a  higher 
sampling  rate.  This  might  reveal  distinguishing  frequency  content  data  at  higher 
frequencies.  If  scintillation  frequencies  can  be  determined,  appropriate  filters  can  be 
constructed  to  detect  scintillation  in  real  time. 

3.  Determine  Solar  Cycle  Dependence 

The  GPS  data  collected  here  were  at  solar  cycle  minimum.  This  is  the  best- 
case  scenario,  as  scintillation  becomes  worse  as  the  solar  cycle  increases  toward 
maximum.  To  determine  how  the  magnitude  of  scintillation  varies  with  changes  in 
the  solar  cycle,  monitoring  low-latitude  ionospheric  scintillation  should  be  repeated 
at  mid-cycle  and  cycle  maximum.  This  could  be  accomplished  using  the  JPL  site  in 
Hawaii  and  a  GPS  receiver  at  Christmas  Island. 

4.  Develop  Daytime  Scintillation  Models 

This  report  was  concerned  primarily  with  nighttime  scintillation  effects  at  low 
latitudes.  In  Chapter  VII,  the  equatorial  anomalies  were  shown  to  contribute 
afternoon  scintillation.  This  same  data  set  can  be  analyzed  to  determine  the  extent 
of  scintillation  magnitude  during  the  daytime. 

5.  Real-Time  Nowcasting 

Real-time  nowcasting  can  be  accomplished  using  GPS  receivers.  The  analysis 
conducted  here  can  be  expanded  to  determine  the  optimum  method  for  attaining 
that  capability.  A  limitation  in  this  data  set  is  that  it  is  for  a  comparatively  benign 
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environment  where  the  magnetic  equator  varies  smoothly  and  constantly  with 
respect  to  the  geographic  equator.  For  better  prediction  of  world-wide  behavior,  this 
experiment  should  be  repeated  in  an  area  of  the  world  where  the  environment  is  not 
so  well  behaved.  One  candidate  area  is  South  America  where  a  few  GPS  receivers 
could  be  emplaced  to  analyze  the  low  latitude  scintillation  there.  This  data  set 
would  help  to  generate  a  real-time  nowcasting  model  for  low-latitude  scintillation. 

6.  Sensor  Augmentation 

It  is  apparent  from  the  TEP  prediction  plots  of  Chapter  VI,  that  not  all  the 
controlling  variables  have  been  identified.  The  connection  between  TEP  and  the 
earth's  magnetic  indices  (i.e.,  the  planetary  magnetic  index,  Kp,  and  the  equatorial 
magnetic  index,  DST)  is  unknown  at  present.  These  inputs  may  improve  TEP 
predictions  and  are  possibilities  for  refinements. 

7.  Correlate  All-Sky  Airglow  and  GPS-TEC 

Scintillation  correlation  between  airglow  intensity  and  TEC  scintillation  at 
zenith  offered  only  limited  analysis  capability.  A  more  meaningful  analysis  could  be 
attained  by  correlating,  with  the  GPS  measurements,  airglow  intensity  values  which 
are  obtained  from  the  same  azimuth  and  elevation  as  the  GPS  satellites. 
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APPENDIX  A. 
GEOGRAPHIC  TO  GEOMAGNETIC  COORDINATES 


In  cartesian  coordinates  the  conversion  can  be  computed  from  geographic  to 
geomagnetic  coordinates  by  (Bhavanani  and  Vancour,  1991) 


mag 


0.320158  -0.928599  -0.187626 
0.945388  0.325947     0 
0.061156  -0.177380  0.982240 


*  r 


geog' 


(63) 


where  r^  is  the  cartesian  coordinate  vector  (x^,  y^,  z^)  for  the  magnetic  latitude 
and  longitude  and  rgeog  is  the  cartesian  coordinate  vector  (x^,  ygeog,  zgeog).  The 
cartesian  geographical  coordinate  vector  can  be  found  by 


geog 


geog 


y 


geog 


geog 


cos(0g)*cos((i>g) 

cos(6p*sin(4>p 

0.9966472  *sin(0g) 


(64) 


where  0g  is  the  geographic  latitude  and  <{>g  is  the  geographic  longitude.  Similarly, 
0m  is  used  for  the  geomagnetic  latitude  and  4>m  for  the  geomagnetic  longitude.  The 
numerically  modified  z  component  in  Eq.  64  accounts  for  the  earth's  eccentricity. 
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APPENDIX  B.   MAGNETIC  DECLINATION  ANGLE 

The  magnetic  declination  angle,  A,  is  the  angle,  with  the  point  of  interest 
(0g,4>g)  on  earth  as  the  apex,  formed  by  the  lines  drawn  from  the  point  of  interest  to 
true  north  and  to  magnetic  north  as  shown  in  Figure  84. 


True  North 


Magnetic 
North 


Figure  84  Magnetic  Declination  Angle. 

After  computing  the  magnetic  coordinates  of  the  point  of  interest,  as  detailed  in 
Appendix  A ,  the  magnetic  declination  angle  can  be  computed.  By  incrementing  the 
point  of  interest  a  minuscule  amount  toward  true  north,  finding  the  magnetic 
coordinates  of  the  new  point,  then  comparing  the  angular  changes,  the  declination 
angle  can  be  computed.  This  is  shown  in  Figure  85. 
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Figure  85  Computing  the  Magnetic  Declination  Angle 

From  the  figure,  the  length  of  y  can  be  seen  to  be  the  difference  between  the  two 
longitudes  multiplied  by  the  cosine  of  the  latitude,  i.e., 


y  =  (lonl  -  lon2)  *  cos(latl) 


(65) 


and  that  x  is  the  difference  of  the  two  latitudes, 


x  =  lat2  -  latl, 


(66) 


With  these  relationships  the  magnetic  declination  angle  can  then  be  found  as 


A  =  tan"1 


y 


(67) 
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APPENDIX  C.  SOLAR  MAGNETIC  ANGLE 

As  the  earth  precesses,  its  axis  tilts  with  respect  to  the  ecliptic  plane  of  the 
sun.  The  earth  tilts  toward  the  sun  in  spring  and  summer  months  and  tilts  away  in 
fall  and  winter.  On  the  days  of  equinox,  the  tilt  is  zero  and  the  sun  is  directly  over 
the  earth's  equator.  The  tilt  maxima  occur  on  the  solstices  with  magnitudes  of 
±23.5°,  with  positive  being  when  the  northern  hemisphere  is  tilted  toward  the  sun. 

The  solar  declination  angle  is  the  angle  formed  (at  the  center  of  the  earth) 
between  the  earth's  and  sun's  ecliptic  planes.  The  declination  angle  limits,  as  with 
the  tilts,  are  ±23.5°. 

The  vernal  equinox,  March  21  of  each  year,  which  is  day  80  on  the  Julian 
calendar  in  non-leap  years,  marks  the  transition  between  a  negative  to  a  positive 
declination  angle  .  Assuming  a  sinusoidal  variation,  an  appropriate  equation  for  the 
solar  declination  angle  is 


6  =  23.5  sin 


2n(JD  -  80)' 


365 


(68) 


/ 


where  JD  is  the  Julian  date  of  the  day  of  interest. 

This  equation  for  the  solar  declination  angle  is  based  upon  the  assumption  of 
sinusoidal  tilt  variation.  For  many  applications,  this  equation  yields  acceptable 
accuracy.  For  additional  accuracy,  the  empirically  based  formulation  given  by 
Bourges  (1985)  can  be  entered  into  a  programmable  device  such  as  a  calculator  or 

143 


PC.  A  MATLAB  program  to  compute  the  solar  declination  angle  in  this  manner  is 
given  in  Appendix  D7. 

The  magnetic  declination  angle  is  the  angle  formed  when  great-circle  lines  are 
drawn  from  a  point-of-interest  to  both  the  geographic  and  magnetic  poles.  A 
MATLAB  program  to  compute  the  magnetic  declination  angle  is  given  in  Appendix 
D5. 

When  the  point-of-interest  is  on  the  equator,  the  magnetic  declination  angle 
indicates  not  only  the  angular  offset  of  the  magnetic  from  the  geographic  poles,  but 
also  the  tilt  of  the  magnetic  equator  with  respect  to  the  geographic  equator.  For 
example,  if  the  magnetic  north  pole  is  offset  10°  eastwardly  from  the  geographic 
north  pole  at  some  point-of-interest  on  the  geographic  equator,  the  magnetic  equator 
will  be  likewise  tilted  10°  from  the  geographic  equator. 

If  an  eastward  magnetic  declination  angle  is  defined  as  positive  (and 
westward  negative),  a  comparison  can  be  made  between  the  solar  and  magnetic 
declination  angles.  When  they  are  the  same,  the  sun  will  set  on  the  magnetic 
equator.  The  solar  magnetic  angle  is  the  difference  between  the  solar  declination  and 
the  magnetic  declination  angles.  With  the  magnetic  declination  angle  given  as  10° 
from  above,  on  approximately  Julian  days  106  and  239  the  solar  magnetic  angle  will 
be  zero. 
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APPENDIX  D.  TEP  PREDICTION  PROGRAMS 

The  prediction  programs  developed  in  Chapter  VI  are  listed  in  this  appendix. 
The  programs  are  written  modularly  so  that  each  module  may  be  easily  converted 
to  a  subroutine  for  the  parent  program  in  the  parent  language.  Each  module  is  listed 
in  a  sub-appendix  starting  with  the  main  program  in  sub-appendix  Dl. 
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Dl.  MAIN 

%  File  Name:  main.m 

% 

%  Purpose:     Serves  as  the  controlling  program  to  call  the  subroutines 

%  and  collects  input  data. 

% 

%  Inputs:        year,  Julian  day,  solar  flux,  frequency  of  transmission 

% 

%  Output:      Predicted  opening  time  of  circuit,  closing  time  of  circuit, 

%  path  loss,  and  probability  of  circuit  existing 

% 

%  Author:       R.  Smith 

%  Mod  Date:  26  June  1996 

%Programs  called:    getsites.m,  gtcirc.m,  flndeq.m,  geo2mag.m,  tepchk.m, 
%  decang.m,  opetime.m,  clostime.m,  prob.m,  ploss.m 

%  Initialize 

dtr  =  pi/ 180; 

%  Get  input  data 

year  =  input('Enter  the  current  year  (e.g.,  1996): '); 

jday  =  input  ('Enter  Julian  Day  (e.g.,  319): ') 

flux  =  input  ('Enter  the  10.7  cm  Solar  Flux:  '); 

freq  =  input  ('Enter  the  frequency  of  interest  (in  MHz): '); 

%  Call  getsites  to  enter  the  latitudes  and  longitudes  of  the 
%  transmit  and  receive  sites. 

getsites 

%  Call  gtcirc  to  find  the  great  circle  azimuths  between  the  sites 
%  and  the  distance  between  them. 

gtcirc 

%  Call  findeq  to  get  the  longitude  of  the  equatorial  crossing  point 
%  of  the  path  between  the  two  sites. 

equatorlong  =  findeq(s2n,n2s,latn,lonn,lats,lons)/dtr; 
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%  Convert  site  geographic  coordinates  to  geomagnetic  coordinates 
%  and  get  magnetic  declination  angle 

[maglatn,maglonn,ma]  =  geo2mag(latn,lonn); 
[maglats,maglons,ma]  =  geo2mag(lats,lons); 
[maglateq,magloneq,ma]  =  geo2mag(0,equator_long); 

%  Check  to  see  if  locations  are  within  TEP  limits. 

tepchk 

%  Get  today's  solar  zenith  angle  and  filtered  solar  magnetic  angle 

[fsma,chi]  =  decang(year,jday,ma); 

%  Compute  circuit  opening  time 

opetime 

%  Compute  circuit  closing  time 

clostime 

%  Compute  probability  of  circuit  actually  opening 

prob 

%  Compute  path  loss 

ploss 
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D2.  GETSITES 

%  File  Name:  getsites.m 

% 

%  Purpose:  Prompts  user  for  the  latitudes  and  longitudes  of  the 

%  two  points  of  interest  (transmit  site  and  receive 

%  site).  A  library  of  named  sites  can  be  augmented  so 

%  that  often  used  sites  can  be  entered  by  name. 

% 

%  Inputs:  There  are  three  ways  to  enter  the  latitudes  and 

%  longitudes:  1)  site  name,  2)  decimal  degrees,  and 

%  3)  degrees  and  minutes. 

% 

%  Author:  R.  Smith 

%  Mod  Date:  24  April  1996 

%  Programs  Called:  hr2deg.m 

%Build  library  of  named  sites 

raro=  [-21.22-159.74]; 
khon=  [21.32-157.85]; 
np=[90  0]; 
sp=[-90  0]; 

%Display  information  concerning  orientation  of  data 

linel  =  ['The  following  inputs  assume  north  and  east  as  positive,']  ; 

line2  =  ['south  and  west  as  negative.'] ; 

disp  (linel) 

disp(line2) 

clear  linel;  clear  line2 

%  Inititialize 

ce  =  'n'; 

%  Input  northernmost  latitude  until  correct  (indicated  by  'y') 

while  ce  ~  =  y 

%  Determine  if  data  is  in  degrees,  minutes  (1),  a  named  site  (2), 
%  or  in  degrees  and  tenths  of  degrees. 
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linel  =  ['If  northernmost  latitude  is  in  DDHHMM  enter  1,']; 

disp(linel) 

clear  linel 

hrmin  =  input('site  name  2;  (otherwise  any  other  key): '); 

if  hrmin  ==1; 

lat  =  input('Enter  latitude  in  DD.HHMM: '); 

Northlatitude  =  lat 

clear  Northjatitude 

latn  =  hr2deg(lat); 

ce  =  inputfls  this  the  correct  latitude?  (y/n):  ','s')  ; 

elseif  hrmin  ==2 

loc  =  input  ('Enter  the  name  of  the  location  of  the  northernmost  site:  ','s') 
dummy  =  eval(loc); 
latn  =  dummy(l); 
lonn  =  dummy  (2); 
ce  =  V; 

else 

latn  =  input  ('Enter  latitude  in  DD.  decimal: '); 

Northlatitude  =  latn 

clear  North_latitude 

ce  =  inputfls  this  the  correct  latitude?  (y/n):  ','s')  ; 
end 

%  End  While 
end 

%  Reinitialize  for  longitude 

ce  =  'n'; 

%  If  a  site  name  was  entered  the  longitude  is  already  recorded 
%  so  that  this  while  block  can  be  skipped 

if  hrmin  =  =  2; 
ce  =  y; 
end 

%Input  northernmost  longitude  until  correct  (indicated  by  'y') 

while  ce  ~  =  'y'; 
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linel  =  ['If  northernmost  longitude  is  in  DDHHMM  type  1,']; 

disp(linel) 

hrmin  =  input(  '(otherwise  any  other  key): '); 

if  hrmin  ==  1; 

Ion  =  inputC'Enter  longitude  in  DD.HHMM: '); 

Northlongitude  =  Ion 

clear  Northlongitude 

lonn  =  hr2deg(lon); 

ce  =  input('Is  this  the  correct  longitude?  (y/n):  ','s')  ; 
else 

lonn  =  input('Enter  longitude  in  DD.decimal: '); 

Northlongitude  =  lonn 

clear  Northlongitude 

ce  =  input('Is  this  the  correct  longitude?  (y/n):  ','s')  ; 
end 

%  End  While 
end 

%  Reinitialize  for  southern  point 

ce  =  'n*; 

while  ce  ~=  "yf; 

linel  =  [If  southernmost  latitude  is  in  DDHHMM  type  1,']; 

disp  (linel) 

hrmin  =  input('site  name  2,  (otherwise  any  other  key): '); 

if  hrmin  ==  1; 

lat  =  input('Enter  latitude  in  DD.HHMM: '); 

Southlatitude  =  lat 

clear  Southjatitude 

lats  =  hr2deg(lat); 

ce  =  input('Is  this  the  correct  latitude?  (y/n):  ','s')  ; 

elseif  hrmin  =  =  2 

loc  =  input('Enter  the  name  of  the  location  of  the  southernmost  site:  ','s'); 
dummy  =  eval(loc); 
lats  =  dummy(l); 
Ions  =  dummy(2); 
ce  =  V; 

151 


else 

end 
end 

lats  =  input  ('Enter  latitude  in 
Southlatitude  =  lats 
clear  South_latitude 
ce  =  input  ('Is  this  the  correct 

DD.decimal:  '); 
latitude?  (y/n):  ','s')  ; 

ce  =  'n'; 

if  hrmin  =  = 
end 

=  2  ce  =  V; 

while 

ce  ~ 

=  y 

linel  =  ['If  southernmost  longitude  is  in  DDHHMM,  type  1,'] 

disp(linel) 

hrmin  =  input  ('(otherwise  any  key): '); 

if  hrmin  ==  1; 

Ion  =  input('Enter  longitude  in  DD.HHMM: '); 

Southlongitude  =  Ion 

clear  Southlongitude 

Ions  =  hr2deg(lon); 

ce  =  inputfls  this  the  correct  longitude?  (y/n):  ','s')  ; 

else 

Ions  =  input('Enter  longitude  in  DD.decimal:  *); 

Southlongitude  =  Ions 

clear  Southlongitude 

ce  =  input('Is  this  the  correct  longitude?  (y/n):  ','s')  ; 
end 
end 

clear  lat;  clear  Ion;  clear  ce;  clear  hrmin 

%  Keep  longitudes  between  0  and  360  degrees 

if  lonn  <  0 

lonn  =  lonn  +  360; 
end 
if  Ions  <  0 

Ions  =  Ions  +  360; 
end 
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D3.   GTCIRC 

%  File  Name:  gtcirc.m 

% 

%  Purpose:  Find  the  distance  between  and  the  great  circle 

%  azimuths  between  two  lat/long  positions  on  earth, 

%  specifically  the  transmit  and  receive  sites. 

% 

%  Inputs:  The  northernmost  latitude  (latn),  northernmost 

%  longitude  (lonn),  the  southenmost  latitude  (lats), 

%  and  the  southernmost  longitude  (Ions).  (All  in  degrees) 

% 

%  Outputs:  Great  circle  azimuths  between  the  two  sites 

%  (X  and  Y  in  radians)  and  the  distance  between  them 

%  (Z  in  radians). 

% 

%  Author:  R.  Smith 

%  Mod  Date:  24  April  1996 

%  Initialize 

dtr  =  pi/180; 

%Prevent  divide  by  zero  in  equations  below 

if  (90  -  latn  <  0.001 ) 
n2s  =  pi/2; 
s2n  =  0; 
dist  =  pi/2  -  lats*dtr; 

elseif  (90  +  lats  <  0.001) 

n2s  =  pi/2; 

s2n  =  0; 

dist  =  pi/2  +  latn*dtr; 
else 

%  Use  spherical  trig  to  determine  the  great  circle  azimuths  from 
%  each  point  to  the  other  (X  and  Y)  and  the  distance  between  them 
%  in  radians  (Z). 

cosdist  =  sin(latn*dtr)*sin(lats*dtr); 

cosdist  =  cosdist  +  cos(latn*dtr)*cos(lats*dtr)*cos((lons  -  lonn)*dtr); 
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cosdist  =  max(cosdist,  -1); 
cosdist  =  min  (cosdist,  1); 

dist  =  acos  (cosdist); 

cosn2s  =  (sin(lats*dtr)  -  sin(latn*dtr)*cosdist)/(cos(latn*dtr)*sin(dist)); 
coss2n  =  (sin(latn*dtr)  -  sin(lats*dtr)*cosdist)/(cos(lats*dtr)*sin(dist)); 

cosn2s  =  max(cosn2s,  -1); 
cosn2s  =  min(cosn2s,  1); 
coss2n  =  max(coss2n,  -1); 
coss2n  =  min(coss2n,  1); 

n2s  =  acos(cosn2s); 
s2n  =  acos(coss2n); 

if  (sin((lons  -  lonn)*dtr)  <=  0) 

n2s  =  2*pi  -  n2s; 
else 

s2n  =  2*pi  -  s2n; 
end 

if  s2n  >  2*pi 

s2n  =  s2n  -  2*pi; 
end 

ifs2n  <  (-2*pi) 

s2n  =  s2n  +  2*pi; 
end 

if  n2s  >  2*pi 

n2s  =  n2s  -  2*pi; 
end 

if  n2s  <  (-2*pi) 

n2s  =  n2s  +  2*pi; 
end 

end 
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D4.   FINDEQ 

%  File  Name:  findeq.m 

% 

%  Purpose:  Given  two  points  north  and  south  of  the  equator 

%  finds  the  point  on  the  equator  which  lies  on 

%  the  great  circle  between  the  two  points. 

% 

%  Inputs:  great  circle  azimuth  from  southern  point  to 

%  northern  point  (X),  great  circle  azimuth  from 

%  northern  point  to  southern  point  (Y),  latitude 

%  and  longitude  of  northern  point  (latn  &  lonn), 

%  latitude  and  longitude  of  southern  point  (lats 

%  Ions) 

% 

%  Output:  longitude  of  equatorial  great  circle  intersection 

%  between  the  two  input  points  (newlon) 

% 

%  Mod  Date:  24  April  19^6 

function  [newlon]  =  findeq(X,Y,latn,lonn,lats,lons) 

%Initialize 

dtr  =  pi/ 180; 
newlat  =  pi/4; 

%For  fewest  iterations  migrate  toward  the  equator  from  the  point 
%which  is  closer  to  the  equator.  Determine  which  point  is  closer 
%and  assign  its  latitude  to  the  variable  lat,  its  longitude  to 
%lon  and  the  great  circle  azimuth  toward  the  other  point  as  hdg. 

if  (abs(latn)  >  =  abs(lats)) 

lat  =  lats*dtr; 

Ion  =  lons*dtr; 

hdg  =  X; 
else 

lat  =  latn*dtr; 

Ion  =  lonn*dtr; 

hdg  =  Y; 
end 

%Perform  iterations  to  tolerance  of  within  0.01  radians  away  from 
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%equator. 

while  (  abs(newlat)  >  0.01) 

%Set  distance  to  move  toward  the  equator  to  an  amount  equal 
%to  the  value  of  the  latitude. 

dist  =  abs(lat); 

%Find  new  latitude  after  moving  distance  dist  toward  equator 

sinNewlat  =  sin  (lat)*cos  (dist)  +  cos(lat)*sin(dist)*cos(hdg); 
newlat  =  asin  (sinNewlat); 

%Find  new  longitude  after  moving  distance  dist  toward  equator 

sinDeltaLon  =  sin(dist)*sin(hdg)/cos(newlat); 
DeltaLon  =  asin  (sinDeltaLon); 
newlon  =  Ion  +  DeltaLon; 

%Maintain  0  <  newlon  <  2  pi 

if  (newlon  <  0) 

newlon  =  newlon  +  2*pi; 
end 

if  (newlon  >  2*pi) 

newlon  =  newlon  -  2* pi; 
end 

%Find  azimuth  from  new  point  toward  equator  along  great  circle 

sinNewHdg  =  cos(lat)*sin(hdg)/cos(newlat); 
newhdg  =  asin(sinNewHdg); 

%Ensure  that  new  azimuth  is  reverses  direction  if  the  equator 
%is  overshot. 

if  (lat  >  newlat) 

if  (DeltaLon  >  0) 

newhdg  =  pi  -  newhdg; 
else 

newhdg  =  -pi  -  newhdg; 
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end 
end 

%keep  0  <  newhdg  <  2  pi 

if  (newhdg  >  2*pi) 

newhdg  =  newhdg  -  2* pi; 
end 

if  (newhdg  <  0) 

newhdg  =  newhdg  +  2*pi; 
end 

%Assign  new  values  for  next  iteration 

lat  =  newlat; 
Ion  =  newlon; 
hdg  =  newhdg; 

end 
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D5.  GE02MAG 


%  File  Name: 

% 

%  Purpose: 

% 

% 

% 

% 

%  Note: 

% 

% 

% 

% 

%  Inputs: 

% 

%  Outputs: 

% 

% 

% 

%  Mod  Date: 

%  Author: 


geo2mag.m 

Given  the  geographic  latitude  and  longitude  of  the 

point  of  interest,  the  geomagnetic  latitude  and 
longitude  are  computed.  The  magnetic  declination 
angle  is  also  computed. 

The  declination  angle  is  the  difference  between 
true  north  and  magnetic  north  at  the  point  of  interest. 
This  value  can  be  found  from  a  navigation  chart 
or  a  geographic  to  geomagnetic  conversion  program. 

geographic  latitude  and  longitude  in  degrees  (lat  and  Ion). 

geomagnetic  latitude  and  longitude  in  degrees 

(latm  and  lonm)  and  the  magnetic  declination  angle 
in  degrees  (magangle). 

26  April  1996 
R.  Smith 


function  [latm,lonm,magangle]  =  geo2mag(lat,lon) 

%Initialize 

dtr  =  pi/180; 
bump  =  0.00001; 

%Convert  to  radians 

latr  =  lat*dtr; 
lonr  =  lon*dtr; 

%Bump  the  latitude  to  make  comparison  from  which  to  determine 
%the  magnetic  declination. 

latl  =  latr  +  bump; 

%Convert  from  spherical  to  cartesian.  The  z  component  accounts 
%for  the  earth  eccentricity. 
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x  =  cos(latr)*cos(lonr); 
y  =  cos(latr)*sin(lonr); 
z  =  0.9966472*sin(latr); 

xl  =  cos(latl)*cos(lonr); 
yl  =  cos(latl)*sin(lonr); 
zl  =  0.9966472*sin(latl); 

%Geographic  to  geomagnetic  conversion  matrix 

m  =  [  0.320158  -0.928599  -0.187626 
0.945388  0.325947  0 
0.061156-0.177380  0.982240]; 

xmag  =  m*  [x;y;z]; 

%Convert  back  to  spherical 

x2  =  xmag  .*  xmag; 

theta  =  acos(xmag(3)/sqrt(x2(l)+x2(2)+x2(3))); 

phi  =  atan2(xmag(2),xmag(l)); 

latm  =  (pi/2  -  theta)/dtr; 
lonm  =  phi/dtr; 

xlmag  =  m  *  [xl;  yl;  zl]; 
x2  =  xmag  .*  xmag; 

theta  =  acos(xlmag(3)/sqrt(x2(l)+x2(2)+x2(3))); 
phi  =  atan2(xlmag(2),xlmag(l)); 

lat2  =  (pi/2  -  theta)/dtr; 
lon2  =  phi/dtr; 

magangle  =  atan(  (lonm-lon2)*cos(latm*dtr)/abs((lat2-latm)))/dtr; 
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D6.  TEPCHK 

%  File  Name:  tepchk.m 

% 

%  Purpose:     Determine  if  the  two  sites  are  within  the  geomagnetic 

%  limit  to  make  TEP  probable 

% 

%  Inputs:       magnetic  latitude  and  longitude  of  transmit  and  receive  sites 

% 

%  Author:       R.  Smith 

%  Mod  Date:  20  June  1996 

%  Initialize 

transflag  =  0; 
latflag  =  0; 
lonflag  =  0; 

%  Check  to  see  that  two  sites  are  in  opposite  hemispheres 

if  (sign(maglatn)  ~  =  sign(maglats)) 

transflag  =  1; 
end 

%  Check  latitudes 

latmin  =  min(abs(maglatn),abs(maglats)); 
latmax  =  max(abs(maglatn),abs(maglats)); 

if  ((latmin  >=  14)  &  (latmax  <=  26)) 

latflag  =  1; 
end 

%  Check  longitudes 

abslon  =  abs(maglonn  -  maglons); 
if  abslon  <=  8 

lonflag  =  1; 
end 

%  Set  TEP  flag  as  1  for  sites  within  limits,  0  otherwise 

if  ((transflag  +  latflag  +  lonflag)  =  =  3) 
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else 


end 


tepflag  =  1; 
tepflag  =  0; 
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D7.  DECANG 

%  File  Nameidecang.m 

% 

%  Purpose:  Computes  the  declination  angle  of  the  sun  at  noon  (GMT) 

%  on  the  specified  day.  The  day  of  the  spring  equinox 

%  is  computed  and  saved  as  nO. 

% 

%  Inputs:  year  (year),  Julian  day  of  year  (n),  and  magnetic  declination 

%  angle  (ma) 

% 

%  Outputs:  solar  declination  angle  in  degrees  (d),  value  of  filtered 

%  solar  magnetic  angle  angles  (dec),  and  day  of  spring  equinox 

(nO) 

% 

%  Author:  R.  Smith 

%  Mod  Date:  24  April  1996 

%  This  algorithm  is  based  on  the  article  by  Bourges  (1985) 

function  [fsma,d]  =  decang(year,jday,ma) 

%Initialize 

w  =  360/365.2422; 
dtr  =  pi/180; 
dec  =  []; 
summag  =  0; 

%Normalize  year  to  reference 

ref  =  year  - 1969; 

%Compute  spring  equinox 

nO  =  78.801  +  (.2422  *  ref)  -  floor(ref/4); 

%Compute  the  solar  declination  angle  for  current  and  previous  149  days 

for  n  =  jday-149:jday 

t  =  n  -  .5  -  nO; 

phi  =  w  *  t  *  dtr; 
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d  =  .3723  +  (23.2567*sin(phi))  ; 

d  =  d  +  .1149*sin(2*phi)  -  .1712*sin(3*phi)  ; 

d  =  d  -  .7580*cos(phi)  +  .3656*cos(2*phi)  +  .0201*cos(3*phi); 

%  Find  solar  magnetic  angle  (SMA)  and  cos(SMA) 

sma  =  d  -  ma; 

cossma  =  cos(sma*dtr); 

%  Integrate  cos(SMA)  with  previous  values 

summag  =  summag  +  cossma; 

end 

%  Normalize  integrated  cos  (SMA)  values  (providing  filtering) 

summag  =  summag/150; 

%  Add  the  value  of  today's  cos  (SMA)  and  the  filtered  previous  values 

fsma  =  cossma  +  summag; 
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D8.   OPETIME 

%  File  Name:  opetime.m 

% 

%  Purpose:     compute  the  predicted  opening  time  of  the  TEP  circuit 

% 

%  Inputs:       solar  flux  (flux),  solar  zenith  angle  (chi) 

% 

%  Output:      Circuit  Opening  Time  (ot) 

% 

%  Author:      R.  Smith 

%  Mod  Date:  26  June  1996 

%  Initialize 

dtr  =  pi/180; 

%  Compute  the  opening  time  derived  from  solar  flux 

otf  =  0.0000181*(flux  .~2)  -  0.0123  *  flux  +  9.46; 

%  Compute  the  opening  time  based  on  today's  solar  zenith  angle 

otchi  =  18.61  -  10.76  *  cos(chi  *  dtr); 

%  Compute  opening  time  using  both  flux  and  chi 

otfx  =  otf  .*  otchi; 

ot  =  0.011  *  (otfx  .^2)  -  1.334  *  otfx  +  61.46 
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D9.  CLOSTIME 

%  File  Name:  clostime.m 

% 

%  Purpose:     compute  the  predicted  closing  time  of  the  TEP  circuit 

% 

%  Inputs:       solar  flux  (flux),  solar  zenith  angle  (chi) 

% 

%  Output:      Circuit  Closing  Time  (ct) 

% 

%  Author:      R.  Smith 

%  Mod  Date:  26  June  1996 


%  Initialize 

dtr  =  pi/ 180; 

num  =  flux  *  sin(chi*dti;  +  10; 

ras  =  num/ 100; 

ct  =  exp(ras)  +  23.4 
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D9.  PROB 

%  File  Name:  prob.m 

% 

%  Purpose:     compute  the  probability  of  TEP  circuit  opening 

% 

%  Inputs:       solar  flux  (flux),  filtered  solar  magnetic  angle  (fsma), 

%  magnetic  declination  angle  (ma) 

% 

%  Output:      Circuit  Opening  probability  (cprob) 

% 

%  Author:      R.  Smith 

%  Mod  Date:  26  June  1996 


%  Normalize  sum  of  cosine  of  solar  magnetic  angle  added  to  the 
%  filtered  value  to  have  maximum  of  one  and  minimum  of  zero 

%  Find  maximum  and  minimum  of  sum  which  is  function  of 
%  the  magnetic  declination  angle 

maxsum  =  -0.00000294673148  *  (ma  ~3)  -  0.00009418067769  *  (ma  ~2); 
maxsum  =  maxsum  +  0.00480080226611  *  ma  +  1.95151868361121; 

minsum  =  -0.00026774397313  *  (ma  ~2)  -  0.00936506722977  *  ma  + 
1.88527023661150; 

%  Find  the  normalized  filtered  cos(SMA) 

diff  =  maxsum  -  minsum; 
drop  =  minsum/diff; 

redsma  =  fsma/diff; 
normsma  =  redsma  -  drop; 
normsma  =  min(normsma,l); 
normsma  =  max(normsma,0); 

%  Compute  probability  envelope  from  solar  flux  and  circuit  probability 

env  =  flux/480  +  .6; 
env  =  min(env,l); 
env  =  max(env,0); 
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cprob  =  env  *  normsma; 
cprob  =  min(cprob,0.99); 
cprob  =  max(cprob,0.01); 

%  Set  probability  to  zero  if  frequency  is  outside  allowable  range 

if  (freq  <  2)  +  (freq  >  200) 

cprob  =  0; 
end 

%  Set  probability  to  zero  if  transmitter  or  receiver  location  is 
%  out  of  range 

if  (tepflag  ==0) 

cprob  =  0; 
end 
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D10.  PLOSS 

%  File  Name:  ploss.m 

% 

%  Purpose:     Compute  the  path  loss  of  the  circuit  from  the  transmitte 

%  to  the  receiver. 

% 

%  Inputs:        magnetic  latitudes  and  longitudes  of  the  transmitter  and  receiver, 

%  the  predicted  opening  and  closing  times  of  the  circuit,  and 

the 

%  transmission  frequency 

% 

%  Ouput:       path  loss  in  dB 

% 

%  Author:      R.  Smith 

%  Mod  Date:  26  Jun  1996 

%  Initialize 

tvec  =  []; 
lossvec  =  []; 

%  Compute  geographical  loss  for  northern  location 

maglat  =  abs(maglatn); 

if  maglat  <=  18 

latfn  =  cos  (2  *pi*  (maglat- 18)/ 16); 
else 

latfh  =  cos(2*pi*(maglat-18)/32); 
end 

%  Compute  geographical  loss  for  southern  location 

maglat  =  abs(maglats); 

if  maglat  <=  18 

latfs  =  cos(2*pi*(maglat-18)/16); 
else 

latfs  =  cos(2*pi*(maglat-18)/32); 
end 

%  Compute  longitudinal  difference  loss 
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delong  =  maglonn  -  maglons; 
lonf  =  cos(2*pi*delong/32); 

%  Compute  total  losses  as  function  of  time  for  predicted  circuit 
%  duration  of  path.  Eleven  points  are  computed  during  this  time 

dur  =  ct  -  ot; 

%  Bump  the  duration  inward  from  the  extremes  to  prevent  log(O) 

bumpdur  =  dur  -  2e-6; 

inc  =  bumpdur/ 10; 

optloss  =  93  +  20  *  loglO(freq); 

fort  =  (ot+le-6):inc:(ct-le-6) 

tlf  =  sin(pi*(t  -  ot)/dur); 

loss  =  optloss  - 10  *  loglO(latfh)  - 10  *  loglO(latfs); 

loss  =  loss  - 10  *  loglO(lonf)  - 10  *  loglO(tlf); 

lossvec  =  [lossvec  loss]; 
tvec  =  [tvec  t]  ; 

end 
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APPENDIX  E.   PROBABILITY  OF  SCINTILLATION 


Two  MATLAB  programs  are  included  in  this  appendix.  The  first,  prob.m, 
predicts  the  probability  of  scintillation  occurrence  at  low  latitudes.  The  second 
program,  solar.m,  is  called  from  within  prob.m. 

El.  PROB.M 

%  File  Name:  prob.m 

% 

%  Purpose:  compute  the  probability  of  equatorial  scintillation 

% 

%  Inputs:  solar  flux  (flux),  filtered  solar  magnetic  angle  (fsma), 

%  Julian  Day  (jday),  magnetic  declination  angle  (ma) 

% 

%  Output:  Scintillation  probability  of  occurrence  (sprob) 

% 

%  Author:  R.  Smith 

%  Mod  Date:  22  Feb  1998 

%  Get  inputs 

jday  =  input('Enter  the  Julian  Day  of  interest  =  =  >  '); 
year  =  input('Now  the  year  for  that  day  =  =  >  '); 
flux  =  input('And  the  solar  flux  for  that  day  =  =  >  '); 
ma  =  input  ('Enter  the  magnetic  declination  angle  =  =  >  '); 
[fsma,dec]  =  decang(year,jday,ma); 

%  Normalize  sum  of  cosine  of  solar  magnetic  angle  added  to  the 
%  filtered  value  to  have  maximum  of  one  and  minimum  of  zero 

%  Find  maximum  and  minimum  of  sum  which  is  function  of 
%  the  magnetic  declination  angle 

maxsum  =  -0.00000294673148  *  (ma  ~3)  -  0.00009418067769  *  (ma  ~2); 
maxsum  =  maxsum  +  0.00480080226611  *  ma  +  1.95151868361121; 

minsum  =  -0.00026774397313  *  (ma  ~2)  -  0.00936506722977  *  ma  + 
1.88527023661150; 

%  Find  the  normalized  filtered  cos(SMA) 
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diff  =  maxsum  -  minsum; 
drop  =  minsum/diff; 

redsma  =  fsma/diff; 
normsma  =  redsma  -  drop; 
normsma  =  min(normsma,l); 
normsma  =  max(normsma,0); 

%  Compute  probability  envelope  from  solar  flux  and  circuit  probability 

env  =  flux/480  +  .6; 
env  =  min(env,l); 
env  =  max(env,0); 

sprob  =  env  *  normsma; 
sprob  =  min(sprob,0.99); 
sprob  =  max(sprob,0.01) 
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E2.  SOLAR.M 

function  [fsma,d]  =  solar (yearjday, ma) 

%  File  Name:  solar. m 

% 

%  Purpose:     Computes  the  declination  angle  of  the  sun  at  noon  (GMT) 

%  on  the  specified  day  and  the  filtered  value  for  the  previous  150 

days. 

%      The  day  of  the  spring  equinox  is  computed  and  saved  as  nO. 

% 

%  Inputs:       year  (year),  Julian  day  of  year  (n),  and  magnetic  declination 

%  angle  (ma) 

% 

%  Outputs:     solar  declination  angle  in  degrees  (d),  value  of  filtered 

%  solar  magnetic  angle  (fsma) 

% 

%  Author:      R.  Smith 

%  Mod  Date:  22  Feb  1998 

%  This  algorithm  is  based  on  the  article  by  Bourges  (1985) 

%Initialize 

w  =  360/365.2422; 
dtr  =  pi/180; 
dec  =  []; 
summag  =  0; 

%Normalize  year  to  reference 

ref  =  year  - 1969; 

%Compute  spring  equinox 

nO  =  78.801  +  (.2422  *  ref)  -  floor(ref/4); 

%Compute  the  solar  declination  angle  for  current  and  previous  149  days 

for  n  =  jday-149:jday 

t  =  n  -  .5  -  nO; 
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phi  =  w  *  t  *  dtr; 

d  =  .3723  +  (23.2567*sin(phi))  ; 

d  =  d  +  .1149*sin(2*phi)  -  .1712*sin(3*phi)  ; 

d  =  d  -  .7580*cos(phi)  +  .3656*cos(2*phi)  +  ,0201*cos(3*phi); 

%  Find  solar  magnetic  angle  (SMA)  and  cos(SMA) 

sma  =  d  -  ma; 

cossma  =  cos(sma*dtr); 

%  Integrate  cos(SMA)  with  previous  values 

summag  =  summag  +  cossma; 

end 

%  Normalize  integrated  cos(SMA)  values  (providing  filtering) 

summag  =  summag/150; 

%  Add  the  value  of  today's  cos(SMA)  and  the  filtered  previous  values 

fsma  =  cossma  +  summag; 
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APPENDIX  F.  SCINTILLATION  MAGNITUDE  PREDICTION 


This  appendix  contains  one  MATLAB  program,  scintdist,  which  predicts  the 
magnitude  of  scintillation  in  TEC  units  at  a  given  time  after  sunset  and  a 
geomagnetic  latitude. 

%  Name:     SCINTDIST.M 

%  Purpose:  Compute  predicted  magnitude  at  given  time  and  place 

% 

%  Inputs:  time  after  sunset  (tas)  and  magnetic  latitude  (maglat) 

% 

%  Output:  Predicted  scintillation  magnitude  (scintmag)  in  TEC  units 

%  Author:    R.  Smith 
%  Last  Mod:  26  Feb  98 

%  Get  time  after  sunset  input 

tas  =  input('Enter  time  (in  hours)  after  sunset  (1-7)  =  =  >  '); 

%  Find  amplitude  as  function  of  time 

if  tas  <  1 

error('Minimum  time  is  one  hour  after  sunset"); 
elseif  tas  >  7; 

error('Maximum  time  is  seven  hours  after  sunset'); 
elseif  tas  <=  5.5 

timemag  =  3.5*(l-exp(-(tas  - l)/.9)); 
else 

timemag  =  3.5*exp(-(t  -  5.5)/.3); 
end 

%  Find  peak  latitude  as  function  of  time 

peaklat  =  12  *  cos(2*pi*(tas  -  2)/17.5); 

%  Get  magnetic  latitude  input 

maglat  =  input  ('Enter  the  magnetic  latitude  of  interest  =  =  >  '); 
maglat  =  abs(maglat); 
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%  Comput  magitude  factor  based  upon  latitude 

if  maglat  >=  peaklat 

latfactor  =  cos  (2 *pi*  (peaklat  -  maglat)/ 12); 
else 

latfactor  =  cos(2*pi*  (peaklat  -  maglat)/32); 
end 
latfactor  =  (latfactor  +  l)/2; 

scintmag  =  timemag  *  latfactor 
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